Method and system for estimating the position of a mobile device

ABSTRACT

A method of estimating the location of a mobile device, comprising the steps of: collecting location information; selecting at least one of a plurality of different location methods to provide a location estimate; and providing a location estimate based on the at least one selected location method.

FIELD OF THE INVENTION

The present invention relates to a location method and in particular butnot exclusively to a method of determining the, location of a mobilestation in a wireless communications network.

BACKGROUND OF THE INVENTION

Wireless cellular communication networks are known. In these networks,the area covered by the network is divided into a number of cells. Eachof the cells has associated with it a base transceiver station. The basetransceiver stations are arranged to communicate with mobile deviceslocated in the cells. The mobile devices can take any suitable form andare typically mobile telephones.

The need for efficient accurate positioning of the mobile telephones hasincreased. In the USA, network operators must be able to provide thelocation of mobile subscribers making emergency calls. Similar proposalsare currently being considered in Europe. Additionally, commercialservices such as tracking services (that is the emergency servicementioned above, people locating, valuable assets location etc),finding/guidance services (proximity services such as yellow pages,direction indicators, point of interest locators etc) and notificationservices (targeted advertisements, traffic alerts, weather alerts,bus/train location, guided tours etc) are currently being proposed.

In the GSM (global system for mobile communications) standard, fourpositioning methods have been included: cell identity and timingadvance, time of arrival, enhanced observed time difference (E-OTD) anda method based on GPS (global positioning system technology either asstandalone GPS or assisted GPS).

The time of arrival method is able to locate handsets with standardsoftware but requires the installation of new network elements such aslocation measurement units at every base station site. Theenhanced—observed time difference method requires the installation oflocation management units at every two to five base transceiver stationsand a software modification in the handset. The assisted GPS methodrequires installation of a GPS receiver and possibly also locationmeasurement units besides the integration of a GPS receiver into ahandset. All of these methods require the introduction of a new networkelement or corresponding functionality responsible for locationcalculation called the serving mobile location centre SMLC.

The timely deployment of location services for all users including thosehaving handsets which do not include the necessary software or hardware,requires that measurements already available in cellular networks shouldbe used. These techniques are important as they allow operators andservice providers to start offering location based services to allcustomers with minimal additional costs whilst waiting for more accurateand sophisticated location technologies to be available. From thetechnical point of view, even when technology such as enhanced-observedtime difference and assisted GPS are fully available, network basedsoftware solutions will still be needed as backup methods when the newstandardised solutions fail or when the requested accuracy can be metwith such a method. Network based software techniques can also be usedas an initial guess for the algorithms used to implement one of thestandard solutions in order to improve the accuracy or speed ofconvergence of those algorithms.

SUMMARY OF THE INVENTION

According to a first aspect of the present invention, there is provideda method of estimating the location of a mobile device, comprising thesteps of: collecting location information; selecting at least one of aplurality of different location methods to provide a location estimate;and providing a location estimate based on the at least one selectedlocation method.

According to a second aspect of the present invention, there is provideda system for estimating the location of a mobile device, comprising:means for collecting location information; means for selecting at leastone of a plurality of different location methods to provide a locationestimate; and means for providing a location estimate based on the atleast one selected location method.

BRIEF DESCRIPTION OF DRAWINGS

For a better understanding of the present invention and as to how thesame may be carried into effect, reference will now be made by way ofexample only to the accompanying drawings in which:

FIG. 1 shows wireless cellular network to which embodiments of thepresent invention can be applied.

FIG. 2 shows a mobile device served by three base transceiver stationsin a wireless cellular network as shown in FIG. 1.

FIG. 3 shows a geometric representation of a base transceiver stationand mobile device in a wireless cellular network as shown in FIG. 1.

FIG. 4 shows an example of antenna gain in a transceiver as shown inFIG. 1.

FIG. 5 shows another example of antenna gain in a transceiver as shownin FIG. 4.

FIG. 6 shows approximated radiation patterns as found in transceivers asshown in FIG. 1.

FIG. 7 shows the Okumura-Hata path loss graphs for modelled wirelesscellular networks as shown in FIG. 1.

FIG. 8 shows a geometric estimate for the confidence region of alocation estimate for a serving cell in a wireless network as shown inFIG. 1.

FIG. 9 shows the confidence region of a location estimate provided by aseries of cell's cell identity.

FIG. 10 shows a geometric representation of the coverage of a cell in awireless cellular network as shown in FIG. 1.

FIG. 11 shows a series of plots of the probability density function of alocation estimate for various values of path loss/attenuation in awireless cellular network as shown in FIG. 1.

FIG. 12 shows a plot of I_(i0) against path loss/attenuation in awireless cellular network as shown in FIG. 1.

FIG. 13 shows an example of the serving area of a cell in a wirelesscellular network as shown in FIG. 1.

FIG. 14 shows the Circular crown confidence region for a locationestimate produced from a constant TA value.

FIG. 15 shows a further geometric view of the confidence region for alocation estimate as shown in FIG. 14.

FIG. 16 shows the geometry used to calculate the confidence region in aCI location estimate.

FIG. 17 shows a flow diagram detailing the steps to provide locationestimates in embodiments of the invention.

DETAILED DESCRIPTION OF EMBODIMENTS OF THE PRESENT INVENTION

Reference is made to FIG. 1 which shows schematically a wirelesstelecommunications network 2 to which embodiments of the presentinvention can be applied. The area covered by the network is divided upinto cells 4. Associated with each cell is a base station 6. In otherembodiments of the invention base stations 6 may be typically referredto as base transceiver stations (BTS). The base stations 6 are arrangedto communicate with user equipment 8 via a wireless connection. The userequipment is typically a mobile device such as a mobile telephone, acomputer, a personal digital assistant (PDA) or the like. In otherembodiments of the invention the user equipment may also be typicallyreferred to as a mobile station (MS).

Embodiments of the present invention are arranged to estimate thelocation of a mobile device (MS) and an associated confidence region.The confidence region is the region in which the mobile device (MS) canbe expected to be located with a given probability. In other words, anestimated position of the mobile device (MS) can be calculated but thatestimate will not be 100% accurate. The confidence region defines anarea in which it is possible to be reasonably certain that the mobiledevice (MS) is located.

Embodiments of the present invention are arranged to combine a number oflocation algorithms which will be described in more detail hereinafter.

The location algorithms used in embodiments of the present invention usecell identity (CI), timing advance (TA) and received signal strength(RX) measurements. Preferred embodiments use three types of algorithms:

-   -   1) algorithms based on the cell identity. These are referred to        as CI algorithms;    -   2) algorithms based on the cell identity and timing advance        information. These are referred to as CI+TA algorithms; and    -   3) algorithms based on the cell identity and received signal        strength measurements. These are referred to as CI+RX        algorithms.

These different algorithms will be described in more detail hereinafter.For convenience these algorithms are divided into two categories. Thefirst includes the CI+RX algorithms and the second includes the othertwo types of algorithm, that is CI and CI+TA algorithms.

It should be appreciated that in some networks, for example the GSMnetwork, the cell identity, timing advance and received signal strengthdata is already available and used for other purposes. This means thatat least some of the algorithms can be implemented without requiring anymodification to the existing handsets.

The location algorithms used in embodiments of the present inventionrequire a certain set of simultaneous equations to be solved. Thesimultaneous equations, depending on the statistical assumptions of themeasurements can be non-linear or linear in the unknowns. If theequations are non-linear an iterative method has to be adopted to find asolution while if the simultaneous equations are linear, a solution inclosed form exists. Closed form algorithms are computationally lighterthan iterative algorithms.

CI+RX Algorithms

The class of algorithms based on the cell identity and received signalstrength will now be described.

Reference is made to FIG. 2 which shows the principles of locationestimates based on received signal strength. Received signal levels aremeasurements of the levels of received signals from the serving basetransceiver station (BTS) 6 (the base transceiver station (BTS) Withwhich the mobile device (MS) 8 is associated) and the six strongestneighbours. Of course other numbers of base transceiver stations can beused. These measurements are performed by the mobile device (MS) 8 andreported to the fixed part of the network when the mobile device (MS) 8is in a dedicated mode. In the idle mode, the mobile device 10 measuresthe received signal level from the best server, that is the basetransceiver station (BTS) on which it is camped and the six strongestneighbours. However, these measurements cannot be reported to the fixedpart of the network as there is no connection with any of the basetransceiver stations.

The level received by a mobile device or more precisely attenuation thatthe received signal has experienced, depends on the reciprocal positionof the mobile device and base transceiver station involved. The receivedsignal level from multiple base transceiver stations can then becombined to estimate the location of the mobile device.

Embodiments of the present invention belong to the class of networkbased software solutions (NBSS). In particular, embodiments of thepresent invention described in this section are based on the use ofsignal level measurements and focus on applications where coverageprediction maps are not available. However, embodiments of the presentinvention can be used where coverage prediction maps are available.Coverage prediction maps are used in some methods where the mobiledevice location is estimated as the location on a coverage map where thevalues of predicted signal strengths best match the levels actuallymeasured.

The measurement/information needed to implement the algorithms describedhereinafter can be separated into network configuration parameters suchas: base station coordinates; sector orientation and base transceiverantenna radiation patterns in the case of sectorised cells; and maximumbase transceiver station downlink transmission power, and measurementssuch as the level of signals received by the mobile device. In additionto that data, knowledge of suitable models to link level measurements tothe reciprocal position of the mobile devices to be located and the basetransceiver stations used in the positioning procedure are required.These will be discussed in more detail below.

The location methods outlined hereinafter estimate the mobile devicecoordinates at the location where a certain function has its minimumvalue. The function is obtained by combining level observations whichare estimates of the attenuation experienced by the signals received bythe mobile device and transmitted by the base transceiver stationsinvolved. Level observations are estimated by subtracting from thelevels received by the mobile device, the contribution of the basetransceiver stations antenna radiation patterns (which is dependent onthe reciprocal angular position between the mobile device and the basetransceiver stations), the path loss (which is dependent on the distancebetween the mobile device and the respective base transceiver station)and other constant factors such as the base transceiver stationstransmission power, cable losses, antenna losses etc.

Clearly, the accuracy of the results provided by the algorithm dependson the accuracy of the information available. The accuracy of the resultcan be improved by including accurate definitions of the antennaradiation patterns and finely tuned propagation models. These models, inprinciple, do not need to be expressed in an analytical form but may beincluded in the location algorithms as look-up tables from which, givena certain mobile device to base transceiver station angle of arrival,the antenna gains can be retrieved. However, this further refinement Canbe omitted from embodiments of the present invention.

It is preferred that the network information such as base transceiverstation coordinates, sector information, antenna radiation patterns andother adjustable parameters should be constantly updated.

In preferred embodiments of the present invention, the models used torepresent signal level measurements may be continuously adaptive tochanges in the radio environment, preferably in an automatic way. Oneway to accomplish this is to adjust the models based on statisticalcomparisons between a certain set of collected measurements and thecorresponding true quantities. This however requires the exact positionof the mobile device to be known in order to calculate the truequantities and to compare them with the measurements. When themeasurement models are derived offline, drive tests can be performedcollecting at the same time level measurements and mobile devicelocation information by means of a GPS (Global Positioning System)receiver co-located with the mobile device. There may be othertechniques which are implemented which provide a more accurate locationtechnique, for example, those using an estimated observe time differenceor arrival technique or an assisted GPS location estimate, that allowthe determination of the measurement models on-line.

The algorithms which are described hereinafter are applications of themaximum likelihood principle to estimate the mobile device coordinatesby processing a set of observations from the level measurements.Different maximum likelihood approaches are set out.

The location algorithms use signal level measurements to estimate themobile device coordinates, x and y, along with the variants of the slowfading σ_(u) ² assumed to be equal for all the base stations involved.

As explained above in embodiments of the present invention described inthis section the input data required by the algorithms includes signalstrength measurements collected by the mobile device (MS), networkparameters, and a suitable path-loss law—in other words a suitable modelfor the wireless transmission characteristics. The list of input datafollows.

Average power of the signals received by the MS from N base transceiverstations (BTS's) measured in decibels, P_(r) ^(i):

  (1)

Coordinates of the N measured BTS's in meters, (x^(i),y^(i)):(x ^(i) ,y ^(i))(i=1, . . . , N);[x ^(i) ]=[y ^(i)=m   (2)

Orientation of the sectors of the N base transceiver stations (BTS's)(in the case of sectorized cells) in radians measured counter-clockwisefrom the x direction, φ_(B) _(i):φ_(B) ^(i)(i=1, . . . ,N) [φ_(B) ^(i)]=rad   (3)

Maximum radiated power from N measured BTS's in decibels:P _(t,max) ^(i) =P _(t) ^(i) +G _(r,max)−Losses^(i) ; [P _(t,max)^(i)]=dB   (4)

The maximum radiated power, P_(t,max) ^(i), represents the maximum powerat the output of the i-th BTS antenna in the direction of maximum gain.It includes the transmitted power, P_(t) ^(i), the maximum gain of theBTS transmit antenna, G_(t,max) ^(i), antennas losses, cable losses,etc.; all measured in dB. To simplify the mathematical model, P_(t,max)^(i) includes also the maximum gain of the MS receive antenna,G_(t,max).

Combined Transmit-Receive antenna pattern for the i-th BTS (i=1 , . . .,N) in decibels:AP _(tr) ^(i)(ψ^(i)(x,y))=AP _(r) ^(i)(ψ^(i)(x,y)−φ_(B) ^(i))+AP_(t)(ψ^(i)(x,y)−φ_(M)); [AP _(tr) ^(i)]=dB   (5)

The combined radiation pattern, AP_(tr) ^(i)(x,y)), represents the gainintroduced by the antenna installed at the BTS site and the antennainstalled in the handset. The gain depends on the reciprocal orientationof such antennas. In the embodiment represented by the above equation,φ_(M) is the orientation of the antenna installed at the MS, φ_(B) ^(i)is the orientation of the antenna installed at the i-th BTS andψ^(i)(x,y) is the angle of arrival of the signal transmitted by the i-thBTS and received by the mobile device (MS), measured in radianscounterclockwise from the x direction (see FIG. 3): $\begin{matrix}{{{{{\psi^{i}\left( {x,y} \right)} = {\tan^{- 1}\frac{y^{i} - y}{x^{i} - x}}};}\left\lbrack \psi^{i} \right\rbrack} = {rad}} & (6)\end{matrix}$

The combined radiation pattern AP_(tr) ^(i), is a function of the MScoordinates (x,y). It includes the radiation pattern of the transmitantenna installed at the i-th BTS, AP_(t) ^(i)(θ), and the radiationpattern of the antenna installed at the mobile device (MS), AP_(r)(θ).In embodiments of the invention the antenna at the MS isomni-directional and the orientation of such antenna is unknown; thusAP_(r)(θ)=0 dB and φ_(M)=0 rad.

The path-loss law for the propagation between MS and the i-th BTSmeasured in decibels:PL ^(i)(d ^(i)(x,y)); [PL ^(i)]=dB   (7)

The path-loss PL^(i)(d^(i)(x,y)) represents the attenuation experiencedby the signal transmitted by the i-th BTS as it propagates further awayfrom the transmit antenna. It is expressed as a function of the distancebetween the MS and the i-th BTS, d^(i), which in turn depends on the MScoordinates:d ^(i)(x,y)={square root}{square root over ( )}(x ^(i) −x)²+(y ^(i) −y)²  (8)

The output of the location algorithms comprises:

An Estimate of the MS location in meters:{circumflex over (x)},ŷ  (9)

An Estimate of the variance of the log-normal slow fading in decibels:{circumflex over (σ)}_(u) ²   (10)

As mentioned previously, the algorithms apply maximum likelihoodestimation principles. They estimate the mobile device (MS) location andvariation of the slow fading by minimising certain scalar functions. Thefunctions are obtained by combining level observations which aremeasurements of the attenuation experienced by the signal received bythe mobile device (MS) with the corresponding expected quantities,determined as a combination of the base transceiver stations antennaradiation patterns (which are dependent on the reciprocal angularposition between the mobile station and respective base transceiverstations) and path loss (which is dependent on the distance between themobile station and the base transceiver stations). By minimising acertain cost function, (a cost function is a generally applied term inmathematics for the description of optimization problems), thealgorithms find the value of the unknown parameters which globallyminimise the difference between the observed attenuation and theexpected attenuation. The cost function is a measure of the expectederrors in the estimation of the location of the mobile device (MS).Therefore by minimising the cost function the expected location error isalso minimised. Embodiments of the invention use the followingalgorithms depending upon assumed characteristics of the wirelessenvironment.

Embodiments of the present invention perform the first algorithm A wherethe wireless environment is assumed to comprise correlated slow fadingcharacteristics and where the slow fading is assumed to have equalvariance statistics. In other words the signal transmitted by each basetransceiver station (BTS) to the mobile device (MS) has a similar butnot identical characteristics. This type of assumption is accurate wherethe transmission paths between base transceiver stations (BTS) andmobile device are similar.

Algorithm A: Maximum Likelihood Estimation with Level ObservationsAssuming Correlated Slow Fading and Equal Variance

1. Calculate the i-th level observation, L^(i), by subtracting from thei-th measured received power, P_(r) ^(i), the maximum power radiated bythe i-th BTS, P_(t,max) ^(i):L ^(i) =P _(r) ^(i) −P _(t,max) ^(i) ; i=1, . . . ,N   (11)

L^(i) is the total attenuation experienced by the signal transmitted bythe i-th BTS while propagating toward the MS. The total attenuationdepends on path-loss, gains introduced by BTS antenna and MS antenna,fluctuations of the radio channel, etc.

2. Stack the level observations from N BTS's in vector L:L=[L¹ 1, . . . ,L^(N)]^(T)   (12)

3. Solve the minimization problem: $\begin{matrix}{\begin{bmatrix}{\hat{\sigma}}_{u}^{2} \\\hat{x} \\\hat{y}\end{bmatrix} = {\arg\quad{\min\limits_{\lbrack\begin{matrix}\sigma_{u}^{2} \\x \\y\end{matrix}\rbrack}\quad{F\left( {x,{y;\sigma_{u}^{2}}} \right)}}}} & (13)\end{matrix}$where the cost function F (x,y; σ_(u) ²) is defined as follows:$\begin{matrix}{{{F\left( {x,{y;\sigma_{u}^{2}}} \right)} = {{\ln\quad\sigma_{u}^{2}} + {\ln{{r_{L}\left( {x,y} \right)}}} + {\frac{1}{\sigma_{u}^{2}}\left( {L - {m_{L}\left( {x,y} \right)}} \right\rbrack^{T}{r_{L}^{- 1}\left( {x,y} \right)}\left( {L - {m_{L}\left( {x,y} \right)}} \right\rbrack}}}{and}} & \quad & (14) \\{{m_{L}\left( {x,y} \right)} = \left\lbrack {{\mu_{L}^{1}\left( {x,y} \right)},\ldots\quad,{\mu_{L}^{N}\left( {x,y} \right)}} \right\rbrack^{T}} & \quad & (15) \\{{\mu_{L}^{i}\left( {x,y} \right)} = {{- {{PL}^{i}\left( {d^{i}\left( {x,y} \right)} \right)}} - {{AP}_{tr}^{i}\left( {\psi^{i}\left( {x,y} \right)} \right)}}} & \quad & (16) \\{\left\lbrack {r_{L}\left( {x,y} \right)} \right\rbrack_{ij} = \left\{ {{\begin{matrix}1 & {i = j} \\\rho_{u}^{i,j} & {i \neq j}\end{matrix}\quad i},{j = 1},\ldots\quad,N} \right.} & \quad & (17)\end{matrix}$

ρ_(u) ^(i,j)(x,y) is the cross-correlation of the slow fading affectingthe signals propagating from BTS^(i) and BTS^(i) toward the MS.

Embodiments of the present invention perform the second algorithm Bwhere the wireless environment is assumed to comprise uncorrelated slowfading characteristics and where the slow fading is assumed to haveequal variance statistics. In other words the signal transmitted by eachbase transceiver station (BTS) to the mobile device (MS) has unrelatedcharacteristics. This type of assumption is accurate where thetransmission paths between base transceiver stations (BTS) and mobiledevice (MS) have no similar components.

Algorithm B: Maximum Likelihood Estimation with Level ObservationsAssuming Uncorrelated Slow Fading and Equal Variance

1. Calculate the i-th level observation by subtracting from the i-thmeasured received power, P_(r) ^(i), the maximum power radiated by thei-th BTS, P_(t,max) ^(i):L ^(i) =P _(r) ^(i) −P _(t,max) ^(i) ; i=1, . . . ,N   (18)

2. Stack level observations from N BTS's in vector L:L=[L¹, . . . ,L^(N)]^(T)   (19)

3. Solve the minimization problem: $\begin{matrix}{\begin{bmatrix}\hat{x} \\\hat{y}\end{bmatrix} = {\arg\quad{\min\limits_{{\lbrack\begin{matrix}x \\y\end{matrix}\rbrack} \in D_{xy}}\quad{F\left( {x,y} \right)}}}} & (20)\end{matrix}$where the cost function F(x,y) is defined as follows: $\begin{matrix}{{F\left( {x,y} \right)} = {\sum\limits_{i = 1}^{N}\left( {L^{i} + {{PL}^{i}\left( {x,y} \right)} + {{AP}_{tr}^{i}\left( {x,y} \right)}} \right)^{2}}} & (21)\end{matrix}$and D_(xy) is the domain of existence of x and y. Several possibledefinitions for D_(xy) are given later in the present application.

4. Calculate {circumflex over (σ)}_(u) ² as{circumflex over (σ)}_(u) ² =F({circumflex over (x)},ŷ)   (22)

Algorithm B differs from Algorithm A in the definition of the costfunction. At the basis of this different definition is a different modelfor the slow fading affecting the signals transmitted by two BTS's. InAlgorithm A the fading is assumed correlated while in Algorithm B isassumed uncorrelated; this results in a simpler definition of the costfunction of Algorithm B.

Embodiments of the present invention perform the third algorithm C wherethe wireless environment is assumed to comprise uncorrelated slow fadingcharacteristics and where the slow fading is assumed to have equalvariance statistics. In other words the signal transmitted by each basetransceiver station (BTS) to the mobile device (MS) has unrelatedcharacteristics. This type of assumption is accurate where thetransmission paths between base transceiver stations (BTS) and mobiledevice (MS) have no similar components.

Algorithm C: Maximum Likelihood Estimation with Level DifferenceObservations Assuming Uncorrelated Slow Fading and Equal Variance.

1. Calculate the i-th level observation by subtracting from the i-thmeasured received power, Pt, the maximum power radiated by the i-th BTS,P_(t,max) ^(i):L ^(i) =P _(r) ^(i) −P _(t,max) ; i=1, . . . ,N   (23)

2. Calculate the j-th level difference observation by subtracting thej-th level observation from the level observation L¹ taken as reference:D ^(i) =L ¹ −L ^(j) ; j=2, . . . ,N   (24)

3. Stack the N−1 difference of level observations in a vector D:D=[D², . . . ,D^(N)]^(T)   (25)

4. Solve the minimization problem $\begin{matrix}{\begin{bmatrix}\hat{x} \\\hat{y}\end{bmatrix} = {\arg\quad{\min\limits_{{\lbrack\begin{matrix}x \\y\end{matrix}\rbrack} \in D_{xy}}\quad{{F\left( {x,y} \right)}\quad{where}}}}} & (26) \\{{F\left( {x,y} \right)} = {{\sum\limits_{j = 2}^{N}\left( {D^{j} - {\mu_{D}^{j}\left( {x,y} \right)}} \right)^{2}} - {\frac{1}{N}\left( {{\sum\limits_{j = 2}^{N}D^{j}} - {\mu_{D}^{j}\left( {x,y} \right)}} \right)^{2}{and}}}} & (27) \\{{\mu_{D}^{j}\left( {x,y} \right)} = {{- \left\lbrack {{{PL}^{1}\left( {d^{1}\left( {x,y} \right)} \right)} - {{PL}^{j}\left( {d^{j}\left( {x,y} \right)} \right)}} \right\rbrack} - \left\lbrack {{{AP}_{tr}^{1}\left( {\psi^{1}\left( {x,y} \right)} \right)} - {{AP}_{tr}^{j}\left( {\psi^{1}\left( {x,y} \right)} \right)}} \right\rbrack}} & (28)\end{matrix}$D_(xy) is the domain of existence of x and y. Several possibledefinitions for selecting D_(xy) are given later.

5. Calculate {circumflex over (σ)}_(u) ² as{circumflex over (σ)}_(u) ² =F({circumflex over (x)},ŷ)   (29)

Algorithm C differs from Algorithms A and B in the definition of theobservations and, as a consequence, in the definition of the costfunction. The observations considered in Algorithm C are differences inthe attenuation experienced by signals transmitted by two differentBTS's and received by the MS. To simplify Algorithm C, moreover, theslow fading processes affecting signals transmitted by two BTS's areassumed uncorrelated, analogously as in Algorithm B.

Algorithms A, B and C described hereinbefore are non-linear algorithms.The above location algorithms are hereafter further detailed anddescribed.

FIG. 3 shows the basic geometry of the problem. FIG. 3 comprises amobile device (MS) 8, and the i-th base transceiver station (BTS) 6. TheMS 8 and BTS 6 exist in a region defined by a Cartesian co-ordinatesystem 301. All angles are defined as being defined from the x-axis inan anti-clockwise direction. The i-th BTS 6 is located at the pointdefined as (x^(i),y^(i)) and arranged to broadcast and receive with amaximum gain direction 303 on an angle defined as φ_(B) ^(i). The mobiledevice MS 8 is located at the point (x,y) and is arranged also totransmit and receive with a maximum gain direction 305 on an angleφ_(M). The MS 8 is positioned in relation to the i-th BTS 6 by a line307 comprising length d¹(x,y) and by an angle ψ^(i)(x,y). In furtherembodiments of the invention alternative co-ordinate systems to theCartesian system are used. In other embodiments of the invention thepolar reference system centred on BTS^(i) with the distance between MSand BTS^(i), d^(i), defined as the radial coordinate and the Angle OfArrival of the signal received by the mobile device, ψ^(i), defined asthe angular coordinate:

The distance between the mobile device (MS) and the i-th basetransceiver station (BTS) is defined according to equation (30):d ^(i)(x,y)={square root}{square root over ((x ^(i) −x)²+(y ^(i)−y)²)}  (30)

Angle of arrival of the signal transmitted by the i-th BTS and receivedby the mobile device (MS) is defined according to equation (31):$\begin{matrix}{{\psi^{i}\left( {x,y} \right)} = {\tan^{- 1}\frac{y^{i} - y}{x^{i} - x}}} & (31)\end{matrix}$

The transformation between x,y (Cartesian) and d^(i), ψ^(i) (radial)coordinates is completed by the following formulas: $\begin{matrix}\left\{ \begin{matrix}{x = {x^{i} + {d^{i}\cos\quad\psi^{i}}}} \\{y = {y^{i} + {d^{i}\sin\quad\psi^{i}}}}\end{matrix} \right. & (32)\end{matrix}$

The location algorithms as described in the above embodiments estimatethe MS location by processing certain level observations. Suchobservations are obtained from signal strength measurements performed bythe MS. This section derives a model for level observations from a modelof signal strength measurements performed by the MS. Such measurementsare, by definition, estimates of the average power of the receivedsignals. A general model to calculate the average power received by themobile at a particular location, (x,y), from the i-th base station isthe following (all quantities are in dB):P _(r) ^(i)(x,y)=P _(t) ^(i) +G _(r) ^(i)(x,y)+G _(t) ^(i)(x,y)−PL^(i)(x,y)−Losses^(i) +u ^(i)(x,y) dB   (33)Where

-   -   P_(r) ^(i)(x,y) is the power of the signal received by the MS at        location (x,y);    -   P_(t) ^(t) is the BTS transmitting power;    -   G_(r) ^(i)(x,y) is the MS antenna gain in the direction of the        i-th BTS;    -   G_(t) ^(i)(x,y) is the antenna gain of the i-th BTS in the        direction of the MS;    -   PL^(i)(x,y) is the path-loss determined by the propagation path        between MS and BTS;    -   the term “Losses^(i)” takes into account the losses due to        antenna feeder, cables, duplex, divider, etc.    -   u^(i)(x,y) is the shadow fading affecting the signal transmitted        by the i-th BTS. It is generally modelled as a random variable        with log-normal distribution (i.e, u^(i)(x,y) measured in dB is        a Gaussian random variable) with standard deviation σ_(u) ^(i).        Typical values for σ_(u) ^(i) range from 5 to 10 decibels:        u ^(i)(x,y)˜N(0,σ_(u) ^(i)) [u ^(i)]=dB   (34)

Equation (33) is the starting point to study techniques for locatingmobiles by means of observations derived from observed level (orreceived level RXLEV) measurements. It is a general enough model toallow introduction of several parameters affecting the propagation suchas antenna radiation patterns, path loss and random fluctuations.

Antenna gains are usually expressed analytically as a function of anangle, θ, which describes the angular distribution of the power radiatedby the antenna when it is connected to a transmitter. In case ofdirectional antennas, the direction of maximum gain is identified by thedirection θ=0:G(θ=0)=G _(max)   (35)

The antenna gain G(θ) can be separated in the sum of a constant termdefining the maximum antenna gain, G_(max), and the radiation pattern,AP(θ), which describes the angular distribution of the radiated power:G(θ)=G _(max) −AP(θ); AP(θ)≧0 dB; AP(θ=0)=0 dB   (35)

FIG. 4 shows a graph of an antenna gain for an antenna as used in anembodiment of the invention. The graph of FIG. 4 comprises the y-axisrepresenting antenna gain G(θ) and the x-axis representing theorientation from the angle of maximum gain. The antenna gain issymmetrical about the θ=0 line. The antenna gain graph comprises a mainlobe (or beam) 401 and four side lobes (or beams) 403, 405, 407, 409either side of the main lobe. The main lobe 401 comprises a maximum gainG_(m) centred at θ=0, The main lobe antenna gain decreases rapidly tozero either side of θ=0. The side lobes 403, 405, 407, 409 comprisesmaller maximum gains, and arranged so that the first side lobe 403 (theside lobe directly adjacent to the main lobe 401) has a larger maximumgain than the second side lobe 405 (the side lobe adjacent to the firstside lobe 403 and the third side lobe 407). The second lobe 405, inturn, has a larger maximum gain than the third lobe 407 (the side lobeadjacent to the second side lobe 405 and the fourth side lobe 409), andthe third side lobe 407 has a larger maximum gain than the fourth sidelobe 409 (the side lobe adjacent to the third side lobe 407).

When analytic expressions of the antennas radiation patterns are notavailable an approximation is needed.

The orientation of the BTS antenna, φ_(B) ^(i) (as shown in FIG. 3)specifies the direction toward which the maximum power is radiated. Thismeans that, given the (x,y) MS coordinates, ψ^(i)(x,y) (e.g., the anglein which BTS′ “sees” the MS) can be calculated and the downlinktransmission gain toward the MS can be expressed as follows:G _(t) ^(i)(x,y)=G _(t) ^(i)(ψ^(i)(x,y)−φ_(B) ^(i))=G_(t,max) ^(i) −AP_(t) ^(i)(ψ^(i)(x,y)−φ_(B) ^(i))   (36)

In principle, if the MS is equipped with a directional antenna and theorientation of the antenna is known, an analogous formula can be used todescribe the gain of the receive antenna at the MS side:G _(r) ^(i)(x,y)=G _(r) ^(i)(ψ^(i)(x,y)−φ_(M))=G _(r,max) −AP_(r)(ψ^(i)(x,y)−φ_(M))   (37)

Although mobile devices can be equipped with directional antennas, it isnot realistic to assume that the orientation of the MS antenna is known;thus in other embodiments of the present invention the MS antennaradiation pattern can be modelled with an averaged omni-directionalpattern by imposing the following constraintsG _(r,max) =G _(t,avg) ; AP _(r)(θ)=0 dB   (38)

The substitution of (36) and (37) in (33) results in the followingcompact expression for the received power:P _(r) ^(i)(x,y)=P _(t,max) ^(i) −AP _(tr) ^(i)(ψ^(i)(x,y))−PL^(i)(x,y)+u ^(i)(x,y)   (39)where the auxiliary definitions P_(t,max) ^(i) and AP_(tr) ^(i) havebeen defined previously in equations (4) and (5) respectively.

Several analytical expressions for the path-loss PL^(i) in equation (33)have been proposed in the scientific literature. In fact, this term isat the basis of the propagation loss prediction models based on whichcellular operators design their networks. In order to calculateprecisely the attenuation experienced by a signal travelling from a basestation to a mobile, finely tuned prediction models including alsoinformation on topography and morphology of the environment should beconsidered. However, when terrain maps are not available, simplifiedmodels must be used. Such available propagation models are known in theart and in embodiments of the present invention the propagation modelthe path loss as a function of the MS-to-BTS distance only with nodependence from the angle of arrival (AOA):PL ^(i)(x,y)=PL ^(i)(d ^(i)(x,y))   (40)

With the above considerations in mind, the average power received by themobile device (MS) at a particular location, (x,y), from the i-th basestation P_(r) ^(i) can be expressed as follows:P _(r) ^(i)(d ^(i)(x,y),ψ^(i)(x,y))=P _(t,max) ^(i) −AP _(tr)^(i)(ψ^(i)(x,y))−PL ^(i)(d ^(i)(x,y))+u ^(i)(x,y)   (41)where previous definitions for; the maximum radiated power from the i-thbase transceiver station P_(t,max) ^(i) [equation (4)], the combinedtransmit-receive antenna pattern AP_(tr) ^(i) [equation (5)], the pathloss between the mobile device (MS) and the i-th base transceiverstation (BTS) PL^(i) [equation (7)] and the log-normal shadow fadingaffecting the signal transmitted by the i-th BTS u^(i) [equation (34)]hold (all quantities are in dB):

Moving from the model as described by equation (41), the hereafterderive and characterize statistically the observations that can bederived from level measurements in a mobile radio network and used forMS location purposes as used in embodiments of the present invention.

One measurement, or observation, as used in the embodiments of theinvention for location purposes in algorithms A, B and C, is thedifference between the average power received by the mobile, P_(r) ^(i),and the maximum power radiated by the i-th BTS, P_(t,max) ^(i):L ^(i) =P _(r) ^(i) −P _(t,max) ^(i)   (42)

This definition is justified by the fact that the MS location-dependentinformation is embedded in the difference P_(r) ^(i)−P_(t,max) ^(i) andthat the maximum radiated power P_(t,max) ^(i) is a parameter reasonablyeasy to calculate, given the network configuration.

According to the propagation model (41), L^(i)(d^(i)(x,y),ψ^(i)(x,y)) isa random variable, due to the stochastic nature of u^(i)(x,y). Sinceu^(i)(x,y) is a Gaussian random variable, L^(i)(d^(i)(x,y),ψ^(i)(x,y))is a Gaussian random variable as well:L^(i)(x,y)˜N(μ_(L) ^(i)(x,y),σ_(L) ^(i)(x,y))   (43)

The mean value μ_(L) ^(i)(x,y) and standard deviation σ_(L) ^(i)(x,y)can be derived as follows: Mean Value of L^(i)(x,y) $\begin{matrix}{{\mu_{L}^{i}\left( {x,y} \right)} = {{E\left\lbrack {L^{i}\left( {x,y} \right)} \right\rbrack} = {{- {{PL}^{i}\left( {d^{i}\left( {x,y} \right)} \right)}} - {{AP}_{tr}^{i}\left( {\psi^{i}\left( {x,y} \right)} \right)}}}} & (44)\end{matrix}$Variance of L^(i)(x,y)(σ_(L) ^(i)(x,y))² =E[(L ^(i)(x,y))²]−(μ_(L) ^(i)(x,y))² =E[(u^(i)(x,y))²]=(σ_(u) ^(i))²   (45)

The mean value and standard deviations of the i-th observation dependson the MS coordinates (x,y). The probability density function (pdf) ofL^(i) conditioned by x and y can be then expressed as $\begin{matrix}{{f_{{L^{i}❘x},y}\left( {{L^{i}❘x},y} \right)} = {\frac{1}{\sqrt{2\pi}{\sigma_{L}^{i} \cdot \left( {x,y} \right)}}\exp\left\{ {- \frac{\left( {L^{i} - {\mu_{L}^{i}\left( {x,y} \right)}} \right)^{2}}{2\left( {\sigma_{L}^{i} \cdot \left( {x,y} \right)} \right)^{2}}} \right\}}} & (46)\end{matrix}$

In embodiments of the present invention for MS location estimationpurposes, level observations from different BTS's are used. For thisreason, the covariance of two level observations is of interest. Thecovariance of two level observation in turn depends on thecross-correlation of the slow fading processes affecting the propagationof a signal from different BTS's. Often the slow fading processes areconsidered uncorrelated. However, as is known a certain correlationexists between signals sent from different BTS's. Therefore inembodiments of the present invention use a model for thecross-correlation of the slow fading, defined as follows:$\begin{matrix}{\frac{E\left\lbrack {{u^{i}\left( {x,y} \right)}{u^{j}\left( {x,y} \right)}} \right\rbrack}{\sigma_{u}^{i}\sigma_{u}^{j}} = {\rho_{u}^{i,j}\left( {x,y} \right)}} & (47)\end{matrix}$

The Covariance of L^(i)(x,y) and L^(j)(x,y) can therefore be shown asthe following (x and y are neglected for notational convenience).E[(L ^(i)−μ_(L) ^(i))(L ^(j)−μ_(L) ^(j))]=E[(L ^(i) L ^(j))]−μ_(L)^(i)μ_(L) ^(j)=σ_(u) ^(i)σ_(u) ^(j)ρ_(u) ^(i,j)   (48)

Level measurements from multiple BTS's, BTS¹, . . . ,BTS^(N) can becollected and used in further embodiments of the invention to estimatethe MS location. Single level observations as shown in equation (42) canbe stacked in a N×1 vector of observations, L:L=[L¹, . . . ,L^(N)]^(T)   (49)which has a multivariate Gaussian distribution:L˜N(m_(L)(x,y),R_(L)(x,y))   (50)

The mean value and covariance matrix of L can in embodiments of thepresent invention be readily calculated by the use of the results fromthe previous equations: The mean value of L(x,y) can be written as:$\begin{matrix}{{m_{L}\left( {x,y} \right)} = \left\lbrack {{\mu_{L}^{1}\left( {x,y} \right)},\ldots\quad,{\mu_{L}^{N}\left( {x,y} \right)}} \right\rbrack^{T}} & (51)\end{matrix}$

The covariance matrix of L(x,y) can be written as:R _(L)(x,y)=E{LL ^(T) }−m _(L) m _(L) ^(T)   (52)

The generic element of R_(L)(x,y) being (see equation (48))$\begin{matrix}{\left\lbrack {R_{L}L\quad\left( {x,y} \right)} \right\rbrack_{ij} = \left\{ \begin{matrix}\left( \sigma_{u}^{i} \right)^{2} & {i = j} \\{\sigma_{u}^{i}\sigma_{u}^{j}{\rho_{u}^{i,j}\left( {x,y} \right)}} & {i \neq j}\end{matrix} \right.} & (53)\end{matrix}$

The probability density function of L conditioned by x and y cantherefore be written as equation (54) (where |R_(L)(x,y)| indicates thedeterminant of the covariance matrix R_(L)(x,y)) $\begin{matrix}{{f_{{L❘x},y}\left( {{L❘x},y} \right)} = {\frac{1}{\left( {2\pi} \right)^{N/2}{{R_{L}\left( {x,y} \right)}}^{1/2}}\exp\left\{ {{- {\frac{1}{2}\left\lbrack {L - {m_{L}\left( {x,y} \right)}} \right\rbrack}^{T}}{{R_{L}^{- 1}\left( {x,y} \right)}\left\lbrack {L - {m_{L}\left( {x,y} \right)}} \right\rbrack}} \right\}}} & (54)\end{matrix}$

In further embodiments of the invention, one simplifying assumption canbe made by considering the slow fading from different BTS's as havingthe same variance: $\begin{matrix}{\sigma_{u}^{i} = {\sigma_{u}^{j}\overset{\Delta}{=}\sigma_{u}}} & (55)\end{matrix}$

This assumption is accurate where the propagation takes place in ahomogeneous communications environment. In other words where thecommunications environment is consistent and similar. Thus the slowfading affecting different BTS's have statistically the same properties.The assumption of the same covariance for all slow fading links resultsin a modified structure of the covariance matrix R_(L)(x,y). Thestructure of the covariance matrix becomes the product between the slowfading variance (which is common to all BTS's), σ_(u) ², and the matrix,r_(L)(x,y), which depends only on the location-dependent crosscorrelations, ρ_(u) ^(i,j)(x,y):σ_(u) ^(i)=σ_(u) ^(j)

σ_(u) R _(L)(x,y)=σ_(u) ² r _(L)(x,y)   (56) $\begin{matrix}{\left\lbrack {r_{L}L\quad\left( {x,y} \right)} \right\rbrack_{ij} = \left\{ \begin{matrix}1 & {i = j} \\{\rho_{u}^{i,j}\left( {x,y} \right)} & {i \neq j}\end{matrix} \right.} & (57)\end{matrix}$

Using the assumption defined above, the probability density function ofL conditioned by x and y in equation (54) can be written as:$\begin{matrix}{{f_{{L❘x},y}\left( {{L❘x},y} \right)} = {\frac{1}{\left( {2\pi} \right)^{N/2}\sigma_{u}{{r_{L}\left( {x,y} \right)}}^{1/2}}\exp\left\{ {{- {\frac{1}{2\sigma_{u}^{2}}\left\lbrack {L - {m_{L}\left( {x,y} \right)}} \right\rbrack}^{T}}{{r_{L}^{- 1}\left( {x,y} \right)}\left\lbrack {L - {m_{L}\left( {x,y} \right)}} \right\rbrack}} \right\}}} & (58)\end{matrix}$

In further embodiments, the slow fading processes are assumed not onlyto have equal variances but to also be uncorrelated. In such anembodiment r_(L)(x,y) becomes an N×N identity matrix, independent of theMS coordinates:r _(L)(x,y)=r _(L) =I   (59)and the probability density function (58) changes in $\begin{matrix}{{f_{{L❘x},y}\left( {{L❘x},y} \right)} = {\frac{1}{\left( {2\pi} \right)^{N/2}\sigma_{u}}\exp\left\{ {{- \frac{1}{2\sigma_{u}^{2}}}{{L - {m_{L}\left( {x,y} \right)}}}^{2}} \right\}}} & (60)\end{matrix}$where (see equations (51) and (44)) $\begin{matrix}{{{L - {m_{L}\left( {x,y} \right)}}}^{2} = {{\left\lbrack {L - {m_{L}\left( {x,y} \right)}} \right\rbrack^{T}\left\lbrack {L - {m_{L}\left( {x,y} \right)}} \right\rbrack}\quad = {\sum\limits_{i = 1}^{N}\left( {L^{i} + {{PL}^{i}\left( {x,y} \right)} + {{AP}_{tr}^{i}\left( {x,y} \right)}} \right)^{2}}}} & (61)\end{matrix}$A Model for the Level Difference Observation

The i-th level observation L defined in equation (42) is obtained fromthe level measured by the MS on the signal transmitted by the i-th BTS.In a further embodiment of the present invention observations can alsobe derived from the difference between the level measured by the mobiledevice from one BTS taken as reference and the level measured by themobile device from another BTS.

Using difference in levels help in eliminating unknown common biases inthe absolute level measurements. Following an approach similar to thosediscussed above, the following shows the statistical characterization ofthe level difference observation.

The difference in level observation is defined in a further embodimentof the present invention as the difference between the level observationfrom the i-th BTS and the level observation from a reference BTS (in thefollowing identified by the index i=1) as defined in (42):D ^(i) =L ¹ −L ^(i) ; i=2, . . . ,N   (62)

D^(i) is a Gaussian variable, due to the stochastic nature of thelog-normal fading:D^(i)(x,y)˜N(μ_(D) ^(i)(x,y),σ_(D) ^(i)(x,y))   (63)

The mean value and covariance of D^(i)(x,y) can be derived as follows:

The mean value of D^(i)(x,y) $\begin{matrix}{{\mu_{D}^{i}\left( {x,y} \right)} = {{E\left\lbrack {D^{i}\left( {x,y} \right)} \right\rbrack} = {{E\left\lbrack {{L^{1}\left( {x,y} \right)} - {L^{i}\left( {x,y} \right)}} \right\rbrack}\quad = {\left( {{\mu_{L}^{1}\left( {x,y} \right)} - {\mu_{L}^{i}\left( {x,y} \right)}} \right)\quad = {{- \left\lbrack {{{PL}^{1}\left( {d^{1}\left( {x,y} \right)} \right)} - {{PL}^{i}\left( {d^{i}\left( {x,y} \right)} \right)}} \right\rbrack} - \quad\left\lbrack {{{AP}_{tr}^{1}\left( {\psi^{1}\left( {x,y} \right)} \right)} - {{AP}_{tr}^{i}\left( {\psi^{i}\left( {x,y} \right)} \right)}} \right\rbrack}}}}} & (64)\end{matrix}$

The covariance of D^(i)(x,y) and D^(j)(x,y)

The covariance between a pair of level difference observations is neededin embodiments of the present invention because the location algorithmsuse joint multiple level difference observations. The general definitionof covariance between D^(i)(x,y) and D^(j)(x,y) is (x and y areneglected to simplify the notation): $\begin{matrix}{{E\left\lbrack {\left( {D^{i} - \mu_{D}^{i}} \right)\left( {D^{j} - \mu_{D}^{j}} \right)} \right\rbrack} = \left\{ \begin{matrix}{\left( \sigma_{D}^{i} \right)^{2};} & {i = j} \\{{{E\left\lbrack {D^{i}D^{j}} \right\rbrack} - {\mu_{D}^{i}\mu_{D}^{j}}};} & {i \neq j}\end{matrix} \right.} & (65)\end{matrix}$

The term resulting when i=j is the variance of the i-th differenceobservation, (σ_(D) ^(i))²: $\begin{matrix}\begin{matrix}{\left( \sigma_{D}^{i} \right)^{2} = {{E\left\lbrack \left( D^{i} \right)^{2} \right\rbrack} - \left( \mu_{D}^{i} \right)^{2}}} \\{= {{E\left\lbrack \left( L^{1} \right)^{2} \right\rbrack} + {E\left\lbrack \left( L^{i} \right)^{2} \right\rbrack} - {2{E\left\lbrack {L^{1}L^{i}} \right\rbrack}} - \left( \mu_{D}^{i} \right)^{2}}} \\{= {\left( \mu_{L}^{1} \right)^{2} + \left( \sigma_{L}^{1} \right)^{2} + \left( \mu_{L}^{i} \right)^{2} + \left( \sigma_{L}^{i} \right)^{2} - {2\sigma_{u}^{1}\sigma_{u}^{i}\rho_{u}^{1,i}} -}} \\{{2\mu_{L}^{1}\mu_{L}^{i}} - \left( {\mu_{L}^{1} - \mu_{L}^{i}} \right)^{2}} \\{= {\left( {\sigma_{u}^{1} - \sigma_{u}^{i}} \right)^{2} + {2\sigma_{u}^{1}{\sigma_{u}^{i}\left( {1 - \rho_{u}^{1,i}} \right)}}}}\end{matrix} & (66)\end{matrix}$where ρ_(u) ^(1,i) is the cross-correlation between the slow fadingaffecting the propagation from BTS¹ and BTS^(i) defined in equation(47).

The term for i≠j can be calculated as follows: $\begin{matrix}\begin{matrix}{{{E\left\lbrack {D^{i}D^{j}} \right\rbrack} - {\mu_{D}^{i}\mu_{D}^{j}}} = {{E\left\lbrack {\left( {L^{1} - L^{i}} \right)\left( {L^{1} - L^{j}} \right)} \right\rbrack} -}} \\{\left( {\mu_{L}^{1} - \mu_{L}^{i}} \right)\left( {\mu_{L}^{1} - \mu_{L}^{j}} \right)} \\{= {{E\left\lbrack \left( L^{1} \right)^{2} \right\rbrack} - {E\left\lbrack {L^{1}L^{i}} \right\rbrack} -}} \\{{E\left\lbrack {L^{1}L^{j}} \right\rbrack} + {E\left\lbrack {L^{i}L^{j}} \right\rbrack} + -} \\{\left( {\mu_{L}^{1} - \mu_{L}^{i}} \right) \cdot \left( {\mu_{L}^{1} - \mu_{L}^{j}} \right)} \\{= {\left( \mu_{L}^{1} \right)^{2} + \left( \sigma_{L}^{1} \right)^{2} - \left( {{\sigma_{u}^{1}\sigma_{u}^{i}\rho_{u}^{1,i}} + {\mu_{L}^{1}\mu_{L}^{i}}} \right) -}} \\{\left( {{\sigma_{u}^{1}\sigma_{u}^{j}\rho_{u}^{1,j}} + {\mu_{L}^{1}\mu_{L}^{j}}} \right) + \left( {{\sigma_{u}^{i}\sigma_{u}^{j}\rho_{u}^{i,j}} + {\mu_{L}^{i}\mu_{L}^{j}}} \right) -} \\{\left( {\mu_{L}^{1} - \mu_{L}^{i}} \right) \cdot \left( {\mu_{L}^{1} - \mu_{L}^{j}} \right)} \\{= {{\sigma_{u}^{1}\left( {\sigma_{u}^{1} - {\sigma_{u}^{i}\rho_{u}^{1,i}} - {\sigma_{u}^{j}\rho_{u}^{1,j}}} \right)} +}} \\{\sigma_{u}^{i}\sigma_{u}^{j}\rho_{u}^{i,j}}\end{matrix} & (67)\end{matrix}$

Thus in embodiments of the present invention the covariance between twolevel difference observations is summarized as: $\begin{matrix}{{E\left\lbrack {\left( {D^{i} - \mu_{D}^{i}} \right)\left( {D^{j} - \mu_{D}^{j}} \right)} \right\rbrack} = \left\{ \begin{matrix}{{\left( {\sigma_{u}^{1} - \sigma_{u}^{i}} \right)^{2} + {2\sigma_{u}^{1}{\sigma_{u}^{i}\left( {1 - \rho_{u}^{1,i}} \right)}}};} & {i = j} \\{{{\sigma_{u}^{1}\left( {\sigma_{u}^{1} - {\sigma_{u}^{i}\rho_{u}^{1,i}} - {\sigma_{u}^{j}\rho_{u}^{1,j}}} \right)} + {\sigma_{u}^{i}\sigma_{u}^{j}\rho_{u}^{i,j}}};} & {i \neq j}\end{matrix} \right.} & (68)\end{matrix}$

The probability density function of D^(i) conditioned by x and y can beexpressed as follows: $\begin{matrix}{{f_{{D^{i}❘x},y}\left( {{D^{i}❘x},y} \right)} = {\frac{1}{\sqrt{2\pi}{\sigma_{D}^{i}\left( {x,y} \right)}}\exp\left\{ {- \frac{\left( {D^{i}{\mu_{D}^{i}\left( {x,y} \right)}} \right)^{2}}{2\left( {\sigma_{D}^{i}\left( {x,y} \right)} \right)^{2}}} \right\}}} & (69)\end{matrix}$

Measurements from multiple BTS's, BTS¹, . . . ,BTS^(N) can be collectedand used to estimate the MS location. In further embodiments of thepresent invention Level difference observations can be stacked in a(N−1)×1 vector of observations, D:D=[D², . . . ,D^(N)]^(T)   (70)which has a multivariate Gaussian distribution:D˜N(m_(D)(x,y),R_(D)(x,y))   (71)

Mean value and covariance matrix of D can be readily calculated from theresults defined above:

The mean value of D(x,y)m _(D)(x,y)=[μ_(D) ²(x,y), . . . ,μ_(D) ^(N)(x,y)]^(T)=[μ_(L)¹(x,y)−μ_(L) ²(x,y), . . . ,μ_(L) ¹(x,y)−μ_(L) ^(N)(x,y)]^(T)   (72)

The covariance matrix of D(x,y)R _(D)(x,y)=E{DD ^(T) }−m _(D) m _(D) ^(T)   (73)

The generic element of R_(D)(X,Y) being (see equation (68)):$\begin{matrix}\begin{matrix}{\left\lbrack {R_{D}\left( {x,y} \right)} \right\rbrack_{ij} = {E\left\lbrack {\left( {D^{i} - \mu_{D}^{i}} \right)\left( {D^{j} - \mu_{D}^{j}} \right)} \right\rbrack}} \\{= \left\{ \begin{matrix}{{\left( {\sigma_{u}^{1} - \sigma_{u}^{i}} \right)^{2} + {2\sigma_{u}^{1}{\sigma_{u}^{i}\left( {1 - \rho_{u}^{1,i}} \right)}}};} & {i = j} \\{{{\sigma_{u}^{1}\left( {\sigma_{u}^{1} - {\sigma_{u}^{i}\rho_{u}^{1,i}} - {\sigma_{u}^{j}\rho_{u}^{1,j}}} \right)} + {\sigma_{u}^{i}\sigma_{u}^{j}\rho_{u}^{i,j}}};} & {i \neq j}\end{matrix} \right.}\end{matrix} & (74)\end{matrix}$

The following probability density function of D conditioned by x and yresults $\begin{matrix}{{f_{{D❘x},y}\left( {{D❘x},y} \right)} = {\frac{1}{\left( {2\pi} \right)^{N/2}{{R_{D}\left( {x,y} \right)}}^{1/2}}\exp\left\{ {{- {\frac{1}{2}\left\lbrack {D - {m_{D}\left( {x,y} \right)}} \right\rbrack}^{T}}{{R_{D}^{- 1}\left( {x,y} \right)}\left\lbrack {D - {m_{D}\left( {x,y} \right)}} \right\rbrack}} \right\}}} & (75)\end{matrix}$

In further embodiments of the present invention the slow fading fromdifferent BTS's are assumed to have the same variance:σ_(u) ^(i)=σ_(u) ^(i)

σ_(u); i,j=2, . . . ,N   (76)

The elements of the covariance matrix R_(D)(X,y) become: $\begin{matrix}{\left\lbrack {R_{D}\left( {x,y} \right)} \right\rbrack_{ij} = \left\{ \begin{matrix}{{2{\sigma_{u}^{2}\left( {1 - \rho_{u}^{1,i}} \right)}};} & {i = j} \\{{\sigma_{u}^{2}\left( {1 - \rho_{u}^{1,i} - \rho_{u}^{1,j} + \rho_{u}^{i,j}} \right)};} & {i \neq j}\end{matrix} \right.} & (77)\end{matrix}$

Moreover, in further embodiments of the invention if the slow fadingprocesses are assumed uncorrelated, in equation (77) ρ_(u) ^(1,i)=ρ_(u)^(1,j) because i,j≠1 and ρ_(u) ^(i,j)=0 when i≠j; thus the genericelement of R_(D) becomes: $\begin{matrix}{\left\lbrack R_{D} \right\rbrack_{ij} = \left\{ \begin{matrix}{{2\sigma_{u}^{2}};} & {i = j} \\{\sigma_{u}^{2};} & {i \neq j}\end{matrix} \right.} & (78)\end{matrix}$

The covariance matrix can be written as $\begin{matrix}{R_{D} = {{\sigma_{u}^{2}\begin{bmatrix}2 & 1 & \cdots & 1 \\1 & 2 & ⋰ & \vdots \\\vdots & ⋰ & ⋰ & 1 \\1 & \cdots & 1 & 2\end{bmatrix}} = {{\sigma_{u}^{2}\left\{ {I + 1} \right\}} = {\sigma_{u}^{2}r_{D}}}}} & (79)\end{matrix}$where I and 1 are two (N−1)×(N−1) matrices: $\begin{matrix}\begin{matrix}{{I = \begin{bmatrix}1 & 0 & \ldots & 0 \\0 & 1 & ⋰ & \vdots \\\vdots & ⋰ & ⋰ & 0 \\0 & \ldots & 0 & 1\end{bmatrix}};} \\{1 = \begin{bmatrix}1 & \ldots & 1 \\\vdots & ⋰ & \vdots \\1 & \ldots & 1\end{bmatrix}}\end{matrix} & (80)\end{matrix}$

Using the assumption of slow fading uncorrelated and with equalvariance, the probability density function (pdf) (75) can be written asfollows $\begin{matrix}{{f_{{D❘x},y}\left( {{D❘x},y} \right)} = {\frac{1}{\left( {2\pi} \right)^{N/2}\sigma_{u}{r_{D}}^{1/2}}\exp\left\{ {{- {\frac{1}{2\sigma_{u}^{2}}\left\lbrack {D - {m_{D}\left( {x,y} \right)}} \right\rbrack}^{T}}{{r_{D}}^{- 1}\left\lbrack {D - {m_{D}\left( {x,y} \right)}} \right\rbrack}} \right\}}} & (81)\end{matrix}$

The determinant of matrix r_(D) and its inverse can be calculatedexplicitly: $\begin{matrix}\begin{matrix}{{{r_{D}} = N};} \\{{r_{D}}^{- 1} = {\frac{1}{r_{D}}\begin{bmatrix}\left( {N - 1} \right) & {- 1} & \ldots & {- 1} \\{- 1} & \left( {N - 1} \right) & ⋰ & \vdots \\\vdots & ⋰ & ⋰ & {- 1} \\{- 1} & \ldots & {- 1} & \left( {N - 1} \right)\end{bmatrix}}} \\{= {I - {\frac{1}{N}1}}}\end{matrix} & (82)\end{matrix}$

Substituting definitions (70) and (72) and the expression of r_(D) ⁻¹just obtained in the argument of the exponent in (81) it results:$\begin{matrix}{{\left\lbrack {D - {m_{D}\left( {x,y} \right)}} \right\rbrack^{T}{{r_{D}}^{- 1}\left\lbrack {D - {m_{D}\left( {x,y} \right)}} \right\rbrack}} = {{\left\lbrack {{D^{2} - \mu_{D}^{2}},\ldots\quad,{D^{N} - \mu_{D}^{N}}} \right\rbrack^{T}{\left\{ {I - \frac{1}{N}} \right\}\begin{bmatrix}{D^{2} - \mu_{D}^{2}} \\\vdots \\{D^{N} - \mu_{D}^{N}}\end{bmatrix}}} = {{\sum\limits_{j = 2}^{N}\quad\left( {D^{j} - {\mu_{D}^{j}\left( {x,y} \right)}} \right)^{2}} - {\frac{1}{N}\left( {{\sum\limits_{j = 2}^{N}\quad D^{j}} - {\mu_{D}^{j}\left( {x,y} \right)}} \right)^{2}}}}} & (83)\end{matrix}$

The location algorithms used in embodiments of the invention areapplications of the Maximum Likelihood (ML) principle to the level andlevel difference observations defined above. The Maximum Likelihoodprinciple as discussed previously is a method widely used in theestimation theory. A brief review of the maximum likelihood principlefollows.

If x εD is an unknown random parameter defined in a certain domain D andy is an observed random parameter. The principle of the MaximumLikelihood (ML) provides an estimate of x by maximizing the jointprobability density function (pdf) of x and y: $\begin{matrix}{{\hat{x}}_{ML} = {\arg{\max\limits_{x \in D}{f\left( {x,y} \right)}}}} & (84)\end{matrix}$

The ML estimates can be also calculated by maximizing the naturallogarithm of f(x,y), Λ(x)=ln f(x,y); usually referred to as thelog-likelihood function: $\begin{matrix}{{\hat{x}}_{ML} = {\arg{\max\limits_{x \in D}{\Lambda(x)}}}} & (85)\end{matrix}$

By writing f(x,y) as the product of the a posteriori probability densityfunction of the observation y given the unknown x, f(y|x), and the priorprobability density function for the unknown x, f(x). In other words theprobability density function of a first unknown x, f(x), multiplied bythe probability density function of the observation y conditional to thefirst unknown x, f(y|x). A similar expansion of the log-likelihoodfunction becomes Λ(x)=ln f(y|x)+ln f(x). In both cases the ML estimatescan be calculated in embodiments of the present invention by solving thefollowing problem: $\begin{matrix}{{\hat{x}}_{ML} = {\arg{\max\limits_{x \in D}\left\{ {{\ln\quad{f\left( {y❘x} \right)}} + {\ln\quad{f(x)}}} \right\}}}} & (86)\end{matrix}$

If no prior probability density function for the unknown x is available,f(x) can be neglected (or equivalently, x can be assumed to be uniformlydistributed over the domain D), resulting in the simplified form ofequation (86) shown in equation (87). $\begin{matrix}{{\hat{x}}_{ML} = {\arg{\max\limits_{x \in D}\left\{ {\ln\quad{f\left( {y❘x} \right)}} \right\}}}} & (87)\end{matrix}$

In the context of MS location with signal level measurements, theobservation y can be either the vector of level observations L definedin (49)y←L=[L ¹ , . . . ,L ^(N)]^(T) ; yεR ^(N×1)   (88)or the vector of level difference observations D defined in (70):y←D=[L ¹ −L ² , . . . ,L ¹ −L ^(N)]^(T) ; yεR ^((N−1)×1)   (89)

The parameter to be estimated x comprises in embodiments of the presentinvention the estimated MS coordinates (x,y) along with the estimatedvariance of the slow fading σ_(u) ²:x=[σ_(u) ²,x,y]^(T)   (90)

It is shown above that both level observations L and level differenceobservations D can be modelled as multivariate Gaussian randomvariables. In such circumstances the ML criterion (84) can be furthersimplified. If the observation y εR^(M×1) conditioned by the unknownparameter x is a M-variate Gaussian random variable with mean valuem_(y) εR^(M×1) and correlation matrix R_(y) εR^(M×M), the probabilitydensity function of observation y conditioned by x is $\begin{matrix}{{f_{y❘x}\left( {y❘x} \right)} = {\frac{1}{\left( {2\pi} \right)^{M/2}{{R_{y}(x)}}^{1/2}}\exp\left\{ {{- {\frac{1}{2}\left\lbrack {y - {m_{y}(x)}} \right\rbrack}^{T}}{{{R_{y}}^{- 1}(x)}\left\lbrack {y - {m_{y}(x)}} \right\rbrack}} \right\}}} & (91)\end{matrix}$where the dependence of R_(y) and m_(y) from x is explicitly indicated.The natural logarithm of f_(y|z)(y|x) is $\begin{matrix}{{\ln\quad{f\left( {y❘x} \right)}} = {{\frac{M}{2}\ln\quad 2\pi} - {\frac{1}{2}\ln{{R_{y}(x)}}} - {{\frac{1}{2}\left\lbrack {y - {m_{y}(x)}} \right\rbrack}^{T}{{{R_{y}}^{- 1}(x)}\left\lbrack {y - {m_{y}(x)}} \right\rbrack}}}} & (92)\end{matrix}$and the ML estimate of x can be calculated according to the ML criterion(89) as follows: $\begin{matrix}{{\hat{x}}_{ML} = {\arg\quad{\min\limits_{x \in \mathcal{D}}\left\{ {{\ln{{R_{y}(x)}}} + {\left\lbrack {y - {m_{y}(x)}} \right\rbrack^{T}{{R_{y}^{- 1}(x)}\left\lbrack {y - {m_{y}(x)}} \right\rbrack}}} \right\}}}} & (93)\end{matrix}$

Equation (93) provides thus the ML criterion for MS location estimationwhen level observations and level difference observations comply withthe Gaussian statistical models derived above. For simplicity, thefollowing example is restricted to the cases in which the variance ofslow fading from different BTS's is assumed equal (the statisticalmodels for the observations are described above).

When N level observations are considered and the variance of the slowfading is assumed equal for all observations, the ML criterion for aGaussian unknown (93) applies with the following definitions for y,m_(y)(x), R_(y)(x) and M in embodiments of the invention:y←L=L¹, . . . ,L^(N)]^(T)m_(y)(x)←m_(L)(x,y)R_(y)(x)←σ_(u) ²r_(L)(x,y)M=N   (94)where${r_{L}\left( {x,y} \right)} = {\frac{1}{\sigma_{u}^{2}}{R_{L}\left( {x,y} \right)}}$and R_(L)(X,y) is the covariance matrix of the level observations L. TheML criterion to estimate x, y and σ_(u) ² is thus: $\begin{matrix}{\begin{bmatrix}{\hat{\sigma}}_{u}^{2} \\\hat{x} \\\hat{y}\end{bmatrix} = {\arg\quad{\min\limits_{\lbrack\begin{matrix}\sigma_{u}^{2} \\x \\y\end{matrix}\rbrack}\left\{ {{\ln\quad\sigma_{u}^{2}} + {\ln{{r_{L}\left( {x,y} \right)}}} + {{\frac{1}{\sigma_{u}^{2}}\left\lbrack {L - {m_{L}\left( {x,y} \right)}} \right\rbrack}^{T}{{r_{L}^{- 1}\left( {x,y} \right)}\left\lbrack {L - {m_{L}\left( {x,y} \right)}} \right\rbrack}}} \right\}}}} & (95)\end{matrix}$

When N level observations are considered and the slow fading is assumeduncorrelated with equal variance for all observations, the ML criterionfor a Gaussian unknown (93) applies with the following definitions fory, m_(y)(x), R_(y)(x) and M: in embodiments of the present invention.y←L=[L¹, . . . ,L^(N)]^(T)m_(y)(x)←m_(L)(x,y)R_(y)(x)←σ_(u) ²IM=N   (96)

The ML criterion to estimate x, y and σ_(u) ² becomes: $\begin{matrix}{\begin{bmatrix}{\hat{\sigma}}_{u}^{2} \\\hat{x} \\\hat{y}\end{bmatrix} = {\arg\quad{\min\limits_{\lbrack\begin{matrix}\sigma_{u}^{2} \\x \\y\end{matrix}\rbrack}\left\{ {{\ln\quad\sigma_{u}^{2}} + {\frac{1}{\sigma_{u}^{2}}{{L - {m_{L}\left( {x,y} \right)}}}^{2}}} \right\}}}} & (97)\end{matrix}$or analogously by inserting the values of ∥L-m_(L)∥ given in equation61. $\begin{matrix}{\begin{bmatrix}{\hat{\sigma}}_{u}^{2} \\\hat{x} \\\hat{y}\end{bmatrix} = {\arg\quad{\min\limits_{\lbrack\begin{matrix}\sigma_{u}^{2} \\x \\y\end{matrix}\rbrack}\left\{ {{\ln\quad\sigma_{u}^{2}} + {\frac{1}{{\sigma\quad}_{u}^{2}}{\sum\limits_{i = 1}^{N}\left( {L^{i} + {{PL}^{i}\left( {x,y} \right)} + {{AP}_{tr}^{i}\left( {x,y} \right)}} \right)^{2}}}} \right\}}}} & (98)\end{matrix}$

In embodiments of the present invention σ_(u) ² can be estimatedseparately. For fixed x and y this estimation of σ_(u) ² is equivalentto finding the (strictly positive) minimum of the function f(s)=lns+K/s,which is s_(min)=K. The estimated value of σ_(u) ² results in$\begin{matrix}{{\hat{\sigma}}_{u}^{2} = {\sum\limits_{i = 1}^{N}\left( {L^{i} + {{PL}^{i}\left( {\hat{x},\hat{y}} \right)} + {{AP}_{tr}^{i}\left( {\hat{x},\hat{y}} \right)}} \right)^{2}}} & (99)\end{matrix}$

The ML estimation of x and y can be found in this embodiment of theinvention by solving the following minimization problem: $\begin{matrix}{\begin{bmatrix}\hat{x} \\\hat{y}\end{bmatrix} = {\arg\quad{\min\limits_{{\lbrack\begin{matrix}x \\y\end{matrix}\rbrack} \in \mathcal{D}_{xy}}\left\{ {\sum\limits_{i = 1}^{N}\left( {L^{i} + {{PL}^{i}\left( {x,y} \right)} + {{AP}_{tr}^{i}\left( {x,y} \right)}} \right)^{2}} \right\}}}} & (100)\end{matrix}$where D_(xy) is the domain of existence of x and y.

The domain D_(xy) can be determined in further embodiments of thepresent invention, for example, by using any of the additional locationinformation available for implementation of Network Based SoftwareSolutions (NBSS's) For instance, if the Cell Identity of the serving BTSis known, D_(xy) can be defined in an embodiment of the invention as thegeographical region served by such BTS. If in addition also the TimingAdvance is available, D_(xy) can be defined in a further embodiment asthe geographical region determined by the intersection of the servingarea of the BTS identified by the given CI. In a further embodimentD_(xy) can be defined as a circular crown or ring with origin at theserving BTS coordinates and inner/outer radii determined according tothe TA value, for example by using the techniques known in the art.Restricting the domain of x and y in embodiments of the invention, forexample according to CI and TA, has two advantages: the first one isthat the TA information is implicitly taken into account; the second oneis that the convergence of the minimization algorithm is made faster.

When N−1 level difference observations are considered and the slowfading is assumed uncorrelated with equal variance for all observations,the ML criterion for Gaussian unknown (93) applies with the followingdefinitions for y, m_(y)(x), R_(y)(x) and M:y←D=[L ¹ −L ² , . . . ,L ¹ −L ^(N)]^(T)m_(y)(x)←m_(D)(x,y)R_(y)(x)←σ_(u) ²r_(D)M=N−1   (101)where r_(D) is the following (N−1)×(N−1) matrix (independent on the MScoordinates): $\begin{matrix}{r_{D} = \begin{bmatrix}2 & 1 & \cdots & 1 \\1 & 2 & ⋰ & \vdots \\\vdots & ⋰ & ⋰ & 1 \\1 & \cdots & 1 & 2\end{bmatrix}_{{({N - 1})} \times {({N - 1})}}} & (102)\end{matrix}$

The ML criterion to estimate x, y and σ_(u) ² can be derived by usingthe result (83): $\begin{matrix}{{\begin{bmatrix}\hat{x} \\\hat{y}\end{bmatrix} = {\arg\quad{\min\limits_{{\lbrack\begin{matrix}x \\y\end{matrix}\rbrack} \in \mathcal{D}_{xy}}\left\{ {{\sum\limits_{j = 2}^{N}\left( {D^{i} - {\mu_{D}^{i}\left( {x,y} \right)}} \right)^{2}} - {\frac{1}{N}{\sum\limits_{j = 2}^{N}\left( {D^{i} - {\mu_{D}^{i}\left( {x,y} \right)}} \right)^{2}}}} \right\}}}}{{\hat{\sigma}}_{u}^{2} = {{\sum\limits_{j = 2}^{N}\left( {D^{j} - {\mu_{D}^{j}\left( {\hat{x},\hat{y}} \right)}} \right)^{2}} - {\frac{1}{N}\left( {{\sum\limits_{j = 2}^{N}D^{j}} - {\mu_{D}^{j}\left( {\hat{x},\hat{y}} \right)}} \right)^{2}}}}{where}} & (103) \\{{D^{j} = {L^{1} - {L^{j}\quad\left( {{j = 2},\ldots\quad,N} \right)}}}{and}} & (104) \\{{\mu_{D}^{j}\left( {x,y} \right)} = {{- \left\lbrack {{{PL}^{1}\left( {d^{1}\left( {x,y} \right)} \right)} - {{PL}^{j}\left( {d^{j}\left( {x,y} \right)} \right)}} \right\rbrack} - \left\lbrack \quad{{{AP}_{tr}^{1}\left( {\psi^{1}\left( {x,y} \right)} \right)} - {{AP}_{tr}^{j}\left( {\psi^{j}\left( {x,y} \right)} \right)}} \right\rbrack}} & (105)\end{matrix}$

D_(xy) is the domain of existence of x and y. Several possibledefinitions for D_(xy) are given later.

Restricting the domain of x and y, for example according to the valuesCI and TA, has two advantages: the first one is that the TA informationis implicitly taken into account in the estimation; the second one isthat the convergence of the minimization algorithm is made faster.

In further embodiments of the present invention alternativeinterpretations of the ML criterion exist. f(x,z) the joint probabilitydensity function of unknown x and observation z can be written asproduct of the conditional probability density function f(x|z) and themarginal probability density function of the observation f(z). Thelog-likelihood function can thus be written as Λ(x)=ln f(x|z)+ln f(z)and the ML estimates is calculated by solving the following problem:$\begin{matrix}{{\hat{x}}_{ML} = {{\arg\quad{\min\limits_{x \in \mathcal{D}}\left\{ {{\ln\quad{f\left( {x❘z} \right)}} + {\ln\quad{f(z)}}} \right\}}} = {\arg\quad{\max\limits_{x \in \mathcal{D}}\left\{ {\ln\quad{f\left( {x❘z} \right)}} \right\}}}}} & (106)\end{matrix}$

The second equality holds because f(z) does not depend on x. If theMaximum Likelihood estimate of x is inside the domain D, then it can becalculated as being the root of the following equation: $\begin{matrix}{{\frac{\partial}{\partial x}\ln\quad{f\left( {x❘z} \right)}} = 0} & (107)\end{matrix}$

The location method used in further embodiments of the present inventionare based on maximum likelihood criterion as defined in (106) where theobservation z comprises CI and TA information from the serving BTS andCI's and RXLEV's (received level) values from all BTS's involved, and xcomprises the unknown coordinates of the MS:x=[x,y]^(T)   (108)

In location service applications, the knowledge of the CI of a certaincell implies that the geographical coordinates of the BTS antenna, aswell as other parameters such as antenna orientation, cell width,transmitted power, etc, are known.

The ML estimate of the MS coordinates is determined in furtherembodiments of the present invention by solving the followingminimization problem: $\begin{matrix}{\left( {\hat{x},\hat{y}} \right) = {\arg\underset{{({x,y})} \in \mathcal{D}}{\quad\max}\quad\ln\quad{f\left( {x,{y❘z}} \right)}}} & (109)\end{matrix}$or, alternatively, by calculating the roots inside the domain D of thefollowing simultaneous equations: $\begin{matrix}\left\{ \begin{matrix}{{\frac{\partial}{\partial x}\ln\quad{f\left( {x,{y❘z}} \right)}} = 0} \\{{\frac{\partial}{\partial y}\ln\quad{f\left( {x,{y❘z}} \right)}} = 0}\end{matrix} \right. & (110)\end{matrix}$

To apply the ML principle, in the form expressed above, the domain ofthe solution D and the conditional probability density function f(x,y|z)need to be determined. These are detailed below.

D is the domain where the location methods described in this documentlook for the solution x[x,y]^(T). D in (109) can be defined by usingsome a priori information on the region where the handset is possiblylocated. Several possibilities exist; four such methods used inembodiments of the invention are described below.

1. D Determined from Cell Identity (CI) of Serving Cell.

In this case, D represents the geographical region where handsetsconnected to the serving cell are most likely located. D can be thusdefined as the confidence region associated to a location estimate basedon the serving cell CI information. A method to determine such aconfidence region is detailed later. FIG. 8 shows such a confidenceregion. The confidence region comprises the difference in area betweentwo, unequal radii, circle segments 801, 803. The two segments 801, 803with a common origin 805 and common arc angles 807 are defined by thefollowing set of parameters: an origin 805 located at point withcoordinates (x_(o),y_(o)), an inner radius R₁, an uncertainty radius R₂,a orientation angle 809 α and inclusion angle 807 β. The orientationangle 809 defines the angle from the x axis to the start of the arc. Theinclusion angle 807 defines the angle from the start of the arc to theend of the arc.

2. D Determined from Cell Identity (CI) and Timing Advance (TA) ofServing Cell.

A method to determine the parameters of the confidence region as shownin FIG. 8 when, both the CI, and the TA information from the serving BTSis available (the TA, in particular, affects the radii R₁ and R₂). Whenthe TA from serving cell is available, D can be thus determined with themethod provided later eventually neglecting the cell sectorization(i.e., assuming α=0 and β=2π).

3. D Determined from Cell Identity (CI) of All Cells Involved inLocation Calculation.

The coordinates of the BTS's involved in the location estimation providethemselves an indication of the geographical region where the MS islocated. D can be thus defined, for instance, by the convex polygonhaving vertices at the coordinates of the outermost BTSs involved in thelocation calculation. The concept is shown in FIG. 9. FIG. 9 shows sixbase transceiver stations (BTS) with the boundary of D defined by fourBTS's 901, 903, 905, 907, and therefore defining a quadrilateral 909 andtwo BTS's 911, 913, located within this D region boundary 909.

4. D Determined from Coverage Prediction Maps.

Following the same criteria proposed in application PCT/EP/01/01147 todetermine the confidence region of a location estimate, coverageprediction maps for serving and/or neighbouring cells can be used todetermine D. If the TA information from the serving BTS is available,the circular crown centred at the serving BTS coordinates with innerradius R₁ and outer radius R₁+R₂ (as shown in FIG. 8) can be used inaddition to coverage maps in determining D. FIG. 8 shows this area asthe area difference between

Expressions for the Probability Density Function f(x,y|z)

Several expressions for the probability density function f(x,y|z) aredetailed below. In embodiments of the present invention the observationz comprises information derived from received level observations(RXLEV).

Given the following definitions:

-   -   N is the number of received level observations (RXLEV's) used to        estimate the MS coordinates;    -   P_(r) ¹, . . . ,P_(r) ^(N) are the received level observations        (RXLEV's) measured by the MS from the N BTS's involved,        expressed in decibel.

P_(t,max) ¹, . . . ,P_(t,max) ^(N) are the maximum radiated power valuesfrom the N BTS's measured in decibels. The i-th maximum radiated power,P_(t,max) ^(i), represents the maximum power at the output of the i-thBTS antenna in the direction of maximum gain. The value P_(t,max) ^(i)as defined in equation (4) comprises transmitted power, P_(t) ^(i),maximum gain of the BTS transmit antenna, G_(t,max) ^(i), antennalosses, cable losses, etc.; all measured in dB. The total (positive)attenuation experienced by the signals transmitted by each of the NBTS's involved while propagating toward the MS can be expressed indecibels asz ^(i) =P _(t,max) ^(i) −P _(r) ^(i)(i=1, . . . ,N)   (111)

The vector z=[z¹, . . . ,z^(N)]^(T) represents the observation based onwhich the MS location is estimated according to criteria (109) or (110).The probability density function of interest here is thusf(x,y|z)=f(x,y|z ¹ , . . . ,z ^(N))   (112)

By using Bayes Theorem, f(x,y|z) can be expressed as follows:$\begin{matrix}\begin{matrix}{{f\left( {x,{y❘z}} \right)} = {{f\left( {{z❘x},y} \right)}\frac{f\left( {x,y} \right)}{f(z)}}} \\{= \frac{{f\left( {{z❘x},y} \right)}{f\left( {x,y} \right)}}{\int_{\infty}^{+ \infty}{\int_{\infty}^{+ \infty}{{f\left( {{z❘x},y} \right)}{f\left( {x,y} \right)}{\mathbb{d}x}{\mathbb{d}y}}}}}\end{matrix} & (113)\end{matrix}$

f(z|x,y) needed in (113) can be evaluated by using the results from theabove section, where the so-called multiple level observation L=[L¹, . .. ,L^(N)]^(T) corresponds to −z defined in the present document. Themost general expression for f(z|x,y)=f_(L|x,y)(−L|x,y), is$\begin{matrix}{{f\left( {{z❘x},y} \right)} = {\frac{1}{\left( {2\pi} \right)^{N/2}{{R_{z}\left( {x,y} \right)}}^{1/2}}\exp\left\{ {{- {\frac{1}{2}\left\lbrack {z - {m_{z}\left( {x,y} \right)}} \right\rbrack}^{T}}{{R_{z}^{- 1}\left( {x,y} \right)}\left\lbrack {z - {m_{z}\left( {x,y} \right)}} \right\rbrack}} \right\}}} & (114)\end{matrix}$where R_(z)(x,y)=E{zz^(T)}−m_(z)m_(z) ^(T)=R_(L)(x,y) is the covariancematrix of z, m_(z)(x,y)=E{z}=[PL¹{d¹(x,y))+AP_(tr) ¹(ψ^(i)(x,y)), . . .,PL^(N)(d^(N)(x,y))+AP_(tr) ^(N)(ψ^(N)(x,y))]^(T)=−m_(L)(x,y) is themean value of z(x,y), PL^(i)(d^(i)(x,y)) and AP_(tr) ^(i)(ψ^(i)(x,y))are the combined Transmit-Receive antenna pattern and path-lossassociated to the signal transmitted by the i-th BTS, respectively,${\psi^{i}\left( {x,y} \right)} = {\tan^{- 1}\frac{y^{1} - y}{x^{i} - x^{=}}}$is the angle of arrival (AOA) of the same signal, d^(i)(x,y)={squareroot}{square root over ((x^(i)−x)²+(y^(i)−y)²)} is the distance betweenthe MS and the i-th BTS and (x¹,y¹), . . . ,(x^(N),y^(N)) are the x,ycoordinates of the N BTS's involved.

In absence of any other a priori information, the MS coordinates can beassumed uniformly distributed in D. This assumption leads to a jointprobability density function of x and y defined as follows:$\begin{matrix}{{f\left( {x,y} \right)} = \left\{ \begin{matrix}\frac{1}{\mathcal{M}(\mathcal{D})} & {\left( {x,y} \right) \in \mathcal{D}} \\0 & {\left( {x,y} \right) \notin \mathcal{D}}\end{matrix} \right.} & (115)\end{matrix}$where M(D) is the size of D, and D can be determined using one of themethods described previously.

Substitution of (115) in (113) leads to the following expression forf(x,y|z): $\begin{matrix}{{f\left( {x,{y❘z}} \right)} = \left\{ \begin{matrix}{\frac{f\left( {{z❘x},y} \right)}{\int{\int_{\mathcal{D}}{{f\left( {{z❘x},y} \right)}{\mathbb{d}x}{\mathbb{d}y}}}}} & {\left( {x,y} \right) \in \mathcal{D}} \\{0} & {\left( {x,y} \right) \notin \mathcal{D}}\end{matrix} \right.} & (116)\end{matrix}$which can be evaluated by using (112).

In a practical implementation it is very difficult to use (116) in theminimization problem (109). For this reason, alternative approximatedefinitions for f(x,y|z) need to be determined. f(x,y|z) can be writtenas product of the probability density functions conditioned by eachsingle measured attenuation f(x,y|z¹), . . . ,f(x,y|z^(N)):$\begin{matrix}{{f\left( {x,{y❘z}} \right)} = {\prod\limits_{i = 1}^{N}{f\left( {x,{y❘z^{i}}} \right)}}} & (117)\end{matrix}$

The i-th probability density function in (117) represents the likelihoodof (x,y) given the attenuation measured from the i-th cell, z^(i).Physically f(x,y|z^(i)) represents the spatial distribution of the MSwhen the signal received by the MS from the i-th BTS experiences anattenuation z^(i). One embodiment of the present invention describedbelow, defines f(x,y|z^(i)) in such a way that the x,y coordinates areuniformly distributed throughout a region, and defined by theobservation z^(i) and/or the radio coverage properties of the i-th cell.A further embodiment of the present invention also described later, doesnot assume the handsets are uniformly distributed but to instead assumesthat the i-th cell is omni-directional, so that the handsets are assumedto be uniformly distributed angularly from the BTS site but notradially. ‘f(x,y z¹): uniform distribution over a region D_(i)

When the MS is assumed to be uniformly distributed in a certaingeographical region D_(i) associated to the i-th cell involved,f(x,y|z¹) has the following expression $\begin{matrix}{{f\left( {x,{y❘z^{i}}} \right)} = \left\{ \begin{matrix}\frac{1}{\mathcal{M}\left( \mathcal{D}_{i} \right)} & {\left( {x,y} \right) \in \mathcal{D}_{i}} \\0 & {\left( {x,y} \right) \notin \mathcal{D}_{i}}\end{matrix} \right.} & (118)\end{matrix}$where M(D_(i)) is the size of the region D_(i). Several possibilitiesexist to determine D_(i);1. D_(i) from Coverage Maps

If coverage maps generated by coverage prediction tools are available,D_(i) in (118) can be defined as the Hearability area of the i-th cell,H_(i). The hearability area identifies the geographical region where thesignals radiated by the i-th BTS reach the handsets with a signalstrength that is above the MS sensitivity level.

The definition (118) takes into account only the identity of the i-thBTS but does not use the actual observed attenuation. One furtherembodiment of the invention refines the definition of f(x,y|z^(i)), byincluding the measured attenuation, is to substitute the hearabilityarea with the coverage area of the i-th cell, N_(i), which identifiesthe geographical region where the signals radiated by the i-th BTS reachthe MS with the attenuation observed from the i-th BTS. From a practicalpoint of view, it could be beneficial to consider a range of values forattenuation instead of a single value; for example z^(i)±Δz^(i), to takeinto account with Δz^(i) of the attenuation's random fluctuations. Inthis case f(x,y|z^(z)) has the same definition as in (118), with H_(i)substituted by N_(i).

2. Analytical Expressions for D_(i)

In absence of coverage maps, D_(i) in (118) may be expressed with manyanalytical functions. One possibility is to define D_(i) as ageneralized sector as shown in FIG. 10. FIG. 10 shows the area definedby D_(i) as the area defined by the area defined by a segment of acircle 1001 with an origin given at coordinates of the i-th BTS(x^(i),y^(i)) 1003, a sector cell orientation 1005, φ^(i), a sectorwidth of twice the angular width 1007 Δφ^(i). and a (front) radius 1011R_(F) ^(i). The area defined is also that included within the area of asmaller circle 1009 with the same origin 1003 and a (back) radius 1013R_(B) ^(i). The area D_(i) can therefore be defined as follows:$\begin{matrix}{\mathcal{D}_{i}\text{:}\quad\left\{ \begin{matrix}{{{\left( {x - x^{i}} \right)^{2} + \left( {y - y^{i}} \right)^{2}} \leq R_{F}^{i}};} & {0 \leq {{{\psi^{i}\left( {x,y} \right)} - \phi^{i}}} \leq {\Delta\phi}^{i}} \\{{{\left( {x - x^{i}} \right)^{2} + \left( {y - y^{i}} \right)^{2}} \leq R_{B}^{i}};} & {{{{\psi^{i}\left( {x,y} \right)} - \phi^{i}}} > {\Delta\phi}^{i}}\end{matrix} \right.} & (119)\end{matrix}$where${\Psi^{i}\left( {x,y} \right)} = {\tan^{- 1}\frac{y^{i} - y}{x^{i} - x}}$is the angle of arrival of the signal transmitted by the i-th BTS.

With the definition (119) for D_(i), M(D_(i)) in (118) is equal to(R_(F) ^(i))²Δφ^(i)+(π−Δφ^(i))(R_(B) ^(i))². In its most generaldefinition, D_(i) represents a sector cell, but can also represent anomni-directional cell, by setting Δφ^(i)=π and R_(B) ^(i)=0.

In an alternate embodiment of the present invention, D_(i) can bedefined as an ellipse with centre at a reference point of coordinates(x_(R) ^(i),y_(R) ^(i)) and with semi-axes σ_(x,i) and σ_(y,i):$\begin{matrix}{{{{\mathcal{D}_{i}\left( {x,y} \right)}\text{:}\quad\frac{\left( {x - x_{R}^{i}} \right)^{2}}{\sigma_{x,i}^{2}}} + \frac{\left( {y - y_{R}^{i}} \right)^{2}}{\sigma_{y,i}^{2}}} \leq 1} & (120)\end{matrix}$

With the definition (120) for D_(i), f(x,y|z^(i)) in (118) is athree-dimensional cylinder with constant height$\frac{1}{M\left( D_{i} \right)} = \frac{1}{{\pi\sigma}_{x,i}\sigma_{y,i}}$and elliptical base. The general elliptical definition allows anapproximation to the coverage of sector cells to be made, but can alsobe applied to omni-directional cells (by setting the origin (x_(R)^(i),y_(R) ^(i)) at the i-th BTS coordinates and σ_(x,i)=σ_(y,i)=R_(F)^(i) in (120), D_(i) is identical to the one obtained with (119) forΔφ^(i)=π and R_(B) ^(i)=0).

-   -   f(x,y|z^(i)): non uniform distribution for omni-directional        cells

In further embodiments of the present invention the probability densityfunction f(x,y|z^(i)) can be calculated from the joint probabilitydensity function of the distance between the MS and the i-th BTS, d^(i),and the angle of arrival from the same BTS ψ^(i)(x,y),f(d^(i),ψ^(i)|z^(i)), as follows: $\begin{matrix}{{f\left( {x,{y❘z^{i}}} \right)} = {{f\left( {{d^{i} = \sqrt{\left( {x - x^{i}} \right)^{2} + \left( {y - y^{i}} \right)^{2}}},{\psi^{i}\quad = {{\tan^{- 1}\frac{y^{i} - y}{x^{i} - x}}❘z^{i}}}} \right)}{{J\left( {x,y} \right)}}}} & (121) \\{{{{{Where}\quad{{J\left( {x,y} \right)}}} = {{❘{\begin{matrix}\frac{\partial d^{i}}{\partial x} & \frac{\partial d^{i}}{\partial y} \\\frac{\partial\psi^{i}}{\partial x} & \frac{\partial\psi^{i}}{\partial y}\end{matrix}❘}}\quad =}}\quad}{\quad\quad\quad}{\quad{{{\quad\quad}❘{\begin{matrix}\frac{x - x^{i}}{\sqrt{\left( {x - x^{i}} \right)^{2} + \left( {y - y^{i}} \right)^{2}}} & \frac{y - y^{i}}{\sqrt{\left( {x - x^{i}} \right)^{2} + \left( {y - y^{i}} \right)^{2}}} \\{- \frac{y - y^{i}}{\sqrt{\left( {x - x^{i}} \right)^{2} + \left( {y - y^{i}} \right)^{2}}}} & \frac{x - x^{i}}{\sqrt{\left( {x - x^{i}} \right)^{2} + \left( {y - y^{i}} \right)^{2}}}\end{matrix}❘}} = \frac{1}{\sqrt{\left( {x - x^{i}} \right)^{2} + \left( {y - y^{i}} \right)^{2}}}}}} & (122)\end{matrix}$

In the following example the case of the omni-directional cell isconsidered. In this example the MS angular coordinate ψ^(i) can beassumed to be independent from the radial coordinate d^(i) and uniformlydistributed about the region [−π,π]. Using these assumptions thefollowing simplified expression for f(d^(i),ψ^(i)|z^(i)) can be defined:$\begin{matrix}{{{f \cdot \left( {d^{i},{\psi^{i}❘z^{i}}} \right)} = {\frac{1}{2\pi}{f \cdot \left( {d^{i}❘z^{i}} \right)}}};{{\left. \psi^{i} \middle| {\leq \pi} \right.,{d^{i} \geq 0}}}} & (123)\end{matrix}$

Strictly speaking, the i-th observed attenuation z^(i) depends ingeneral on the distance between MS and the i-th BTS, but also on thegain of the transmit antenna in the direction of the MS. If the cellsare omni-directional, as assumed here, the antenna gain contribution canbe neglected, since BTS antenna radiate uniformly in all directions andthe attenuation can be considered as the function of the MS-BTS distanceonly. In this sense, the i-th attenuation z is equal to the i-thso-called path-loss PL^(i)(d^(i)), dependent only on the distancebetween MS and the i-th BTS:z ^(i) =PL ^(i)(d ^(i))   (124)

Thus in summary, the probability density function f(x,y|z^(i)), withomni-directional cells can be obtained by inserting (122), (123), and(124) in (121), as follows: $\begin{matrix}{{f\left( {x,{y❘z^{i}}} \right)} = {\frac{1}{2\pi}{f\left( {d^{i} = {\sqrt{\left( {x - x^{i}} \right)^{2} + \left( {y - y^{i}} \right)^{2}}❘{PL}^{i}}} \right)}\frac{1}{\sqrt{\left( {x - x^{i}} \right)^{2} + \left( {y - y^{i}} \right)^{2}}}}} & (125)\end{matrix}$

It has been shown in experimental and theoretical studies that theaverage received signal power decreases logarithmically with thedistance between transmitter and receiver, both in indoor and in outdoorenvironments; thus the path-loss at a distance d≧d_(o) in decibels canbe expressed as $\begin{matrix}{{{{PL}(d)} = {{{PL}\left( d_{0} \right)} + {10n\quad{\log_{10}\left( \frac{\mathbb{d}}{\mathbb{d}_{0}} \right)}} + u}};{d \geq d_{0}}} & (126)\end{matrix}$

where n is the environment-dependent propagation exponent, d_(o) is the“close-in reference distance” and PL(d_(o)) is the average path-lossexperienced at a distance d_(o) from the transmitter. In free space thevalue of n is 2 but the value n grows when the density of obstructionsincreases. Table 1 lists typical path-loss exponents in differentenvironments. TABLE 1 Environment Propagation exponent, n Free space 2Urban area cellular radio 2.7 ÷ 3.5 Shadowed urban area cellular radio 3÷ 5 In-building line of sight 1.6 ÷ 1.8 Obstructed in-building 4 ÷ 6Obstructed in-factories 2 ÷ 3d_(o) must also be selected on the basis of the environment. When cellsare large, d_(o) is usually set to 1 km; in case of micro-cells, the“close in reference distance” is usually smaller (it can vary from 1 mto 100 m). PL(d_(o)) is therefore calculated in embodiments of theinvention using experimental data. When this is not possible, PL(d_(o))can be estimated in further embodiments by using the free-space pathloss law; if d_(o) is close enough to the transmitter the idealizedcondition of propagation in free-space can be assumed (λ=c/f is thesignal wavelength, c is the speed of light and f the frequency):${{PL}_{{free}\quad{space}}\left( d_{0} \right)} = {10\quad{{\log\left( \frac{4\pi\quad d_{0}}{\lambda} \right)}^{2}.}}$

In the model (126) u represents the shadow fading affecting the signaltransmitted by the i-th cell. It is generally modelled as a randomvariable with log-normal distribution with standard deviation au (i.e.,u measured in dB is a Gaussian random variable: u˜N(0, σ_(u))). Typicalvalues for σ_(u) range from 5 to 10 decibel. By definingA=PL(d ₀)−10n log₁₀ d ₀   (127)B=10nthe model (126) can be re-written as.PL(d)=A+B log₁₀ d+u   (128)(A well known model for the path-loss compliant with (128) is theOkumura-Hata model, where d is the distance between MS and BTS measuredin kilometers, A and B are functions of signal frequency, f, basestation effective antenna height, h_(B), mobile terminal antenna height,h_(m), and City Type (either “Large City” or “Small/Medium City”). TheOkumura-Hata model is coherent with the formula (126) provided that d₀=1km, A=PL(d₀), and B=10n.)

For a given path-loss PL, the distance d≧d₀ has the following expression$\begin{matrix}{d = 10^{- \frac{{PL} - A - u}{B}}} & (129)\end{matrix}$

For a given path-loss value and using the assumption of log-normal slowfading, the probability density function of d in (129) can be calculatedwith known standard random variable transformation techniques:$\begin{matrix}{{{{f\left( d \middle| {PL} \right)} = {\frac{B/{C\left( d_{0} \right)}}{\sqrt{2\pi}\sigma_{u}\quad\ln\quad 10}\frac{1}{d}\exp\left\{ {{- \frac{1}{2\sigma_{u}^{2}}}\left( {{B\quad\log_{10}d} - {PL} + A} \right)^{2}} \right\}}};}{d \geq d_{0}}} & (130)\end{matrix}$where C(d_(o)) is a normalization factor introduced so that ∫_(d) ₀^(∞)f(p|PL)dp=1:C(d ₀)=∫_(d) ₀ ^(∞) f(ρ|PL)dp   (131)

Inserting (130) in (125), the probability density function of x and yfor a given observed attenuation from the i-th omni-directional cell(see (124)) is as follows: $\begin{matrix}{{f\left( {x,\left. y \middle| z^{i} \right.} \right)} = {\frac{B^{i}/{C^{i}\left( d_{0} \right)}}{\left( {2\pi} \right)^{3/2}\sigma_{u}^{i}\quad\ln\quad 10}\frac{\exp\left\{ {{- \frac{1}{2\sigma_{u}^{i^{2}}}}\left( {{B^{i}\log_{10}\sqrt{\left( {x - x^{i}} \right)^{2} + \left( {y - y^{i}} \right)^{2}}} - z^{i} + A^{i}} \right)^{2}} \right\}}{\left( {x - x^{i}} \right)^{2} + \left( {y - y^{i}} \right)^{2}}}} & (132)\end{matrix}$valid for d^(i)(x,y)={square root}{square root over((x−x^(i))²+(y−y^(i))²)}≧d₀.

FIG. 11 shows a series of plots of the probability density functionsf(x,y|z^(i)) using equation (132) for different values of attenuationz^(i)=PL^(i). FIG. 11 shows that the probability density function hascircular symmetry around the vertical axis going through the i-th BTScoordinates. It can be seen that, as the path-loss increases, theprobability density function spreads. At low values of attenuation theprobability density function curves peak close to the i-th BTS andrapidly fall towards zero as the distance from the i-th BTS increases.At higher values of attenuation the peaks of the probability densityfunction curves move further away from the i-th BTS coordinates. As thevalue of attenuation increases the distribution becomes flatter—i.e. thepeak probability density function value is smaller but the rate ofincrease and decrease of the probability density function is lower. Thisspreading produces the suggestion that at higher levels of attenuationthe probability of handsets is further away from the BTS site grows.

-   -   f(x,y|z^(i)): empirical Gaussian distribution for        omni-directional cells

In embodiments of the present invention the probability density functionf(x,y|z^(i)) (132) is approximated by a bi-variate Gaussian probabilitydensity function defined as follows: $\begin{matrix}{{{f_{\mathcal{G}}\left( {x,\left. y \middle| z^{i} \right.} \right)} = {\frac{1}{2\pi{R}^{1/2}}\exp\left\{ {{- {\frac{1}{2}\left\lbrack {x - m} \right\rbrack}^{T}}{R^{- 1}\left\lbrack {x - m} \right\rbrack}} \right\}}},} & (133) \\{where} & \quad \\{x = \begin{bmatrix}x \\y\end{bmatrix}} & (134) \\{m = {{E\left\{ x \middle| z^{i} \right\}} = {\begin{bmatrix}{E\left\{ x \middle| z^{i} \right\}} \\{E\left\{ y \middle| z^{i} \right\}}\end{bmatrix} = \begin{bmatrix}\mu_{x,i} \\\mu_{y,i}\end{bmatrix}}}} & (135) \\{R = {{{E\left\{ {xx}^{T} \right\}} - {mm}^{T}}\quad = \begin{bmatrix}{{E\left\{ x^{2} \middle| z^{i} \right\}} - \left( \mu_{x,i} \right)^{2}} & {{E\left\{ {xy} \middle| z^{i} \right\}} - {\mu_{x,i}\mu_{y,i}}} \\{{E\left\{ {xy} \middle| z^{i} \right\}} - {\mu_{x,i}\mu_{y,i}}} & {{E\left\{ y^{2} \middle| z^{i} \right\}} - \left( \mu_{y,i} \right)^{2}}\end{bmatrix}}} & (136)\end{matrix}$

To derive the Gaussian probability density function in (133), m and R(more precisely, its determinant |R| and its inverse R⁻¹) must bedetermined. It can be noticed that the average values of x and y withprobability density function f(x,y|z^(i)) defined in (132) are$\begin{matrix}{{\mu_{x,i} = {{E\left\{ x \middle| z^{i} \right\}} = x^{i}}};{\mu_{y,i} = {{E\left\{ y \middle| z^{i} \right\}} = {y^{i}\quad{thus}}}}} & (137) \\{m = \begin{bmatrix}x^{i} \\y^{i}\end{bmatrix}} & (138)\end{matrix}$

To obtain an expression for the correlation matrix R in (136), theexpected values E{x²|z^(i)}, E{y²|z^(i)}, and E{xy|z^(i)} need to becalculated.

The analytic expression of E{x²|z^(i)} can be obtained as (the secondequality in the following is obtained by solving the integral in polarcoordinates: $\begin{matrix}{{\phi = {\tan^{- 1}\frac{y - y^{i}}{x - x^{i}}}};{\rho = {{\sqrt{\left( {x - x^{i}} \right)^{2} + \left( {y - y^{i}} \right)^{2}}\text{)}}:}}} & (139) \\{{E\left\{ x^{2} \middle| z^{i} \right\}} = {\int{\int_{\sqrt{{({x - x^{i}})}^{2} + {({y - y^{i}})}^{2}} \geq d_{0}}{x^{2}{f\left( {x,\left. y \middle| z^{i} \right.} \right)}{\mathbb{d}x}{\mathbb{d}y}}}}} & \quad \\{\quad{= {\frac{B^{i}/{C^{i}\left( d_{0} \right)}}{\left( {2\pi} \right)^{3/2}\sigma_{u}^{i}\ln\quad 10}{\int_{- \pi}^{+ \pi}{{\mathbb{d}\quad\phi}{\int_{d_{0}}^{+ \infty}{{\mathbb{d}\rho}\quad{\rho\left( {x^{i} + {\rho\quad\cos\quad\phi}} \right)}^{2}}}}}}}} & \quad \\{\quad\frac{\exp\left\{ {{- \frac{1}{2\sigma_{u}^{i^{2}}}}\left( {{B^{i}\log_{10}\rho} - z^{i} + A^{i}} \right)^{2}} \right\}}{\rho^{2}}} & \quad \\{\quad{= {\frac{B^{i}/{C^{i}\left( d_{0} \right)}}{\left( {2\pi} \right)^{3/2}\sigma_{u}^{i}\ln\quad 10}\left( x^{i} \right)^{2}\underset{\underset{2\pi}{︸}}{\int_{- \pi}^{+ \pi}{\mathbb{d}\quad\phi}}{\int_{d_{0}}^{+ \infty}\frac{1}{\rho}}}}\quad} & \quad \\{\quad{{\exp\left\{ {{- \frac{1}{2\sigma_{u}^{i^{2}}}}\left( {{B^{i}\log_{10}\rho} - z^{i} + A^{i}} \right)^{2}} \right\}{\mathbb{d}\rho}} +}} & \quad \\{\quad{\frac{B^{i}/{C^{i}\left( d_{0} \right)}}{\left( {2\pi} \right)^{3/2}\sigma_{u}^{i}\ln\quad 10}2x^{i}\underset{\underset{0}{︸}}{\int_{- \pi}^{+ \pi}{\cos\quad\phi\quad{\mathbb{d}\phi}}}\int_{d_{0}}^{+ \infty}}} & \quad \\{\quad{{\exp\left\{ {{- \frac{1}{2\sigma_{u}^{i^{2}}}}\left( {{B^{i}\log_{10}\rho} - z^{i} + A^{i}} \right)^{2}} \right\}{\mathbb{d}\rho}} +}} & \quad \\{\quad{\frac{B^{i}/{C^{i}\left( d_{0} \right)}}{\left( {2\pi} \right)^{3/2}\sigma_{u}^{i}\ln\quad 10}\underset{\underset{\pi}{︸}}{\int_{- \pi}^{+ \pi}{\cos^{2}\phi\quad{\mathbb{d}\phi}}}{\int_{d_{0}}^{+ \infty}\rho}}} & \quad \\{\quad{\exp\left\{ {{- \frac{1}{2\sigma_{u}^{i^{2}}}}\left( {{B^{i}\log_{10}\rho} - z^{i} + A^{i}} \right)^{2}} \right\}{\mathbb{d}\rho}}} & \quad \\{\quad{= {{\left( x^{i} \right)^{2}{\overset{\sim}{\mathcal{I}}}_{i}} + \mathcal{I}_{i}}}} & \quad\end{matrix}$where Ĩ_(i) and I_(i) represent, respectively, the first and thirdintegrals above: $\begin{matrix}{{{\overset{\sim}{\mathcal{I}}}_{i} = {\frac{1}{2{C^{i}\left( d_{0} \right)}}{{erfc}\left( \frac{{B^{i}\log_{10}d_{0}} - z^{i} + A^{i}}{\sqrt{2}\sigma_{u}^{i}} \right)}}}{\mathcal{I}_{i} = {\frac{1}{4{C^{i}\left( d_{0} \right)}}\exp\left\{ {2\frac{\ln\quad 10}{B^{i}}\left( {z^{i} - A^{i} + \frac{\sigma_{u}^{i^{2}}\ln\quad 10}{B^{i}}} \right)} \right\}{{erfc}\left( \frac{{B^{i}\log_{10}d_{0}} - z^{i} + A^{i} - {2\frac{\sigma_{u}^{i^{2}}\ln\quad 10}{B^{i}}}}{\sqrt{2}\sigma_{u}^{i}} \right)}}}} & (140)\end{matrix}$

The following definition for the complementary error function is used.$\begin{matrix}{{\frac{1}{2}{{erfc}\left( \frac{u}{\sqrt{2}} \right)}} = {\int_{u}^{\infty}{\frac{{\mathbb{e}}^{{- \theta^{2}}/2}}{\sqrt{2\pi}}{\mathbb{d}\theta}}}} & (141)\end{matrix}$

With a very similar derivation, it can be proved that E{y²|z^(i)} hasthe following expression:E{y ² |z ^(i)}=(y ^(i))² {tilde over (L)} _(i) +I _(i)   (142)

Finally, the analytic expression of the correlation term E{xy|z^(i)} canbe obtained: $\begin{matrix}\begin{matrix}{{E\left\{ {xy} \middle| z^{i} \right\}} = {\int{\int_{\sqrt{{({x - x^{i}})}^{2} + {({y - y^{i}})}^{2}} \geq d_{0}}{{{xyf}\left( {x,\left. y \middle| z^{i} \right.} \right)}{\mathbb{d}x}{\mathbb{d}y}}}}} \\{= {\frac{B^{i}/{C^{i}\left( d_{0} \right)}}{\left( {2\pi} \right)^{3/2}\sigma_{u}^{i}\ln\quad 10}{\int_{- \pi}^{+ \pi}{{\mathbb{d}\phi}{\int_{d_{0}}^{+ \infty}{{\mathbb{d}\rho}\quad\rho}}}}}} \\{\left( {x^{i} + {\rho\quad\cos\quad\phi}} \right)\left( {y^{i} + {\rho\quad\sin\quad\phi}} \right)} \\{\frac{\exp\left\{ {{- \frac{1}{2\sigma_{u}^{i^{2}}}}\left( {{B^{i}\log_{10}\rho} - z^{i} + A^{i}} \right)^{2}} \right\}}{\rho^{2}}} \\{= {x^{i}y^{i}\mathcal{I}_{i}}}\end{matrix} & (143)\end{matrix}$

Using the results (139), (142), and (143) in the definition for R givenin equation (136), the following R results $\begin{matrix}{R = {\begin{bmatrix}{\mathcal{I}_{i} + {\left( x^{i} \right)^{2}\left( {{\overset{\sim}{\mathcal{I}}}_{i} - 1} \right)}} & {x^{i}{y^{i}\left( {{\overset{\sim}{\mathcal{I}}}_{i} - 1} \right)}} \\{x^{i}{y^{i}\left( {{\overset{\sim}{\mathcal{I}}}_{i} - 1} \right)}} & {\mathcal{I}_{i} + {\left( y^{i} \right)^{2}\left( {{\overset{\sim}{\mathcal{I}}}_{i} - 1} \right)}}\end{bmatrix}.}} & (144)\end{matrix}$

The determinant of R being|R|=(I_(i))²+(x^(i))²I_(i)(Ĩ_(i)−1)+(y^(i))²I_(i)(Ĩ_(i)−1) and itsinverse being $\begin{matrix}{R^{- 1} = {{\frac{1}{R}\begin{bmatrix}{\mathcal{I}_{i} + {\left( x^{i} \right)^{2}\left( {{\overset{\sim}{\mathcal{I}}}_{i} - 1} \right)}} & {{- x^{i}}{y^{i}\left( {{\overset{\sim}{\mathcal{I}}}_{i} - 1} \right)}} \\{{- x^{i}}{y^{i}\left( {{\overset{\sim}{\mathcal{I}}}_{i} - 1} \right)}} & {\mathcal{I}_{i} + {\left( y^{i} \right)^{2}\left( {{\overset{\sim}{\mathcal{I}}}_{i} - 1} \right)}}\end{bmatrix}}.}} & (145)\end{matrix}$

With the above results it is possible to define the Gaussian probabilitydensity function as shown in equation (139) as follows: $\begin{matrix}\begin{matrix}{{f_{\mathcal{G}}\left( {x,\left. y \middle| z^{i} \right.} \right)} = {\frac{1}{2\pi{R}^{1/2}}\exp\left\{ {\frac{- 1}{2{R}}\left\lbrack {{I_{i}\left( {x - x^{i}} \right)}^{2} + {I_{i}\left( {y - y^{i}} \right)}^{2} +} \right.} \right.}} \\{\left( {{\overset{\sim}{I}}_{i} - 1} \right)\left( {{x^{i}\left( {x - x^{i}} \right)} - {y^{i}\left( {y - y^{i}} \right)}} \right)^{2}}\end{matrix} & (146)\end{matrix}$

In further embodiments of the present invention the use of the path-lossmodel (128) is extended to distances below the close-in distance {i.e.,for d<d_(o)) or, alternatively, the close-in distance is made tend tod_(o)→0 the above functions remain well behaving and the probabilitydensity function tends to zero as d tends to zero. In the limit d_(o)→0,the term Ĩ_(i) tends to 1 (in fact, (0.5)erf{c(−∞)}=1), C^(i)(d₀)→1, and$\begin{matrix}{\left. {R}\rightarrow\left( I_{i0} \right)^{2} \right.;\left. R^{- 1}\rightarrow{{\frac{1}{x_{i0}}\begin{bmatrix}1 & 0 \\0 & 1\end{bmatrix}}\quad\left( d_{0}\rightarrow 0 \right)\quad{where}} \right.} & (147) \\{I_{i0} = {{\lim\limits_{d_{0}\rightarrow 0}I_{i}} = {\frac{1}{2}\exp\left\{ {2\frac{\ln\quad 10}{B^{i}}\left( {z^{i} - A^{i} + \frac{\sigma_{u}^{i^{2}}\ln\quad 10}{B^{i}}} \right)} \right\}}}} & (148)\end{matrix}$

Inserting the above result for d₀→0 into (133), the followingapproximated Gaussian probability density function results:$\begin{matrix}{{f_{\mathcal{G}0}\left( {x,\left. y \middle| z^{i} \right.} \right)}\frac{1}{2\pi\quad I_{i0}}\exp\left\{ {- \frac{\left( {x - x^{i}} \right)^{2} + \left( {y - y^{i}} \right)^{2}}{2I_{i0}}} \right\}} & (149)\end{matrix}$

In FIG. 12, the behaviour of I_(i0) as a function of thepath-loss/attenuation z^(i) is shown. FIG. 12 shows a graph of I_(i0) ofrange 0 to 8 against path loss/attenuation z^(i) of range 110 to 145 dB.The graph plot resembles an exponential type plot with a I_(i0) ofslightly above 0 at z^(i)=110 dB rising initially slowly but increasingrapidly as Z^(i) passes 140 dB.

-   -   f(x,y|z^(i)): non uniform distribution for omni-directional        cells

The location algorithms in further embodiments may be applied to nonuniform probability density functions in omni-directional cells, Bytaking logarithms of the terms from equation (117), the minimizationproblem can be rewritten as $\begin{matrix}{\left( {\hat{x},\hat{y}} \right) = {\arg\quad{\max\limits_{{({x,y})} \in \mathcal{D}}{\sum\limits_{i = 1}^{N}{\ln\quad{f\left( {x,\left. y \middle| z^{i} \right.} \right)}}}}}} & (150)\end{matrix}$or, alternatively, $\begin{matrix}\left\{ \begin{matrix}{{\sum\limits_{i = 1}^{N}{\frac{\partial}{\partial x}\ln\quad{f\left( {x,\left. y \middle| z^{i} \right.} \right)}}} = 0} \\{{\sum\limits_{i = 1}^{N}{\frac{\partial}{\partial y}\ln\quad{f\left( {x,\left. y \middle| z^{i} \right.} \right)}}} = 0}\end{matrix} \right. & (151)\end{matrix}$

As described above, equation (138), defines the probability densityfunction f(x,y|z^(i)) in the environment created by omni-directionalcells, based on the logarithmic path-loss model (128). By using thisdefinition for the probability density function f(x,y|z^(i)), thefollowing expressions for the partial derivatives can be found$\begin{matrix}\left\{ {\begin{matrix}{{\frac{\partial}{\partial x}\ln\quad{f\left( {x,\left. y \middle| z^{i} \right.} \right)}} = {{F^{i}\left( {x,y} \right)}\left( {x - x^{i}} \right)}} \\{{\frac{\partial}{\partial y}\ln\quad{f\left( {x,\left. y \middle| z^{i} \right.} \right)}} = {{F^{i}\left( {x,y} \right)}\left( {y - y^{i}} \right)}}\end{matrix}\quad{where}} \right. & (152) \\{{F^{i}\left( {x,y} \right)} = \frac{2{B^{i}/{C^{i}\left( d_{0} \right)}}}{\left( {2\pi} \right)^{3/2}\sigma_{u}^{i}\ln\quad 10}} & (153) \\{\quad{\frac{\exp\left\{ {{- \frac{1}{2\sigma_{u}^{i^{2}}}}\left( {{B^{i}\log_{10}{d^{i}\left( {x,y} \right)}} - z^{i} + A^{i}} \right)^{2}} \right\}}{\left\lbrack {d^{i}\left( {x,y} \right)} \right\rbrack^{4}} \cdot}} & \quad \\{\quad\left\lbrack \frac{B^{i}\left( {{B^{i}\log_{10}{d^{i}\left( {x,y} \right)}} - z^{i} + A^{i}} \right)}{2\sigma_{u}^{i^{2}}\ln\quad 10} \right\rbrack} & \quad\end{matrix}$

The MS location can be thus estimated in further embodiments of thepresent invention by solving iteratively the following set ofsimultaneous nonlinear equations $\begin{matrix}\left\{ {\begin{matrix}{{\sum\limits_{i = 1}^{N}\quad{{F^{i}\left( {x,y} \right)}\left( {x - x^{i}} \right)}} = 0} \\{{\sum\limits_{i = 1}^{N}{{F^{i}\left( {x,y} \right)}\left( {y - y^{i}} \right)}} = 0}\end{matrix};{\left( {x,y} \right) \in \mathcal{D}}} \right. & (154)\end{matrix}$

Equation 154 defines algorithm D which is a non linear algorithm.

-   -   f(x,y|z^(i)): empirical Gaussian distribution for        omni-directional cells

Earlier two Gaussian approximations to the probability density functionf(x,y|z^(i)) were defined. The first approximation provided wasf_(G)(x,y|z^(i)), in equation (146). By using such expression, thefollowing partial derivatives result: $\begin{matrix}\left\{ {\begin{matrix}{{\frac{\partial}{\partial x}\ln\quad{f\left( {x,\left. y \middle| z^{i} \right.} \right)}} = \left\lbrack {{{- \frac{I_{i}}{R}}\left( {x - x^{i}} \right)} - {\frac{\left( {{\overset{\sim}{I}}_{i} - 1} \right)}{R}\left\{ {{\left( x^{i} \right)^{2}x} - {x^{i}{y^{i}\left( {y - y^{i}} \right)}}} \right\}}} \right\rbrack} \\{{\frac{\partial}{\partial y}\ln\quad{f\left( {x,\left. y \middle| z^{i} \right.} \right)}} = \left\lbrack {{{- \frac{I_{i}}{R}}\left( {y - y^{i}} \right)} - {\frac{\left( {{\overset{\sim}{I}}_{i} - 1} \right)}{R}\left\{ {{\left( y^{i} \right)^{2}y} - {x^{i}{y^{i}\left( {x - x^{i}} \right)}}} \right\}}} \right\rbrack}\end{matrix},} \right. & (155)\end{matrix}$

Where the determinant of R is given by|R|=(I_(i))²+(x^(i))²I_(i)(Ĩ₁−1)′(y^(i))²I_(i)(Ĩ_(i)−1). The MS locationcan be thus estimated in embodiments of the invention by solvingiteratively the following set of simultaneous nonlinear equations:$\begin{matrix}\left\{ {\begin{matrix}{{\sum\limits_{i = 1}^{N}\left\lbrack {{{- \frac{I_{i}}{R}}\left( {x - x^{i}} \right)} - {\frac{\left( {{\overset{\sim}{I}}_{i} - 1} \right)}{R}\left\{ {{\left( x^{i} \right)^{2}x} - {x^{i}{y^{i}\left( {y - y^{i}} \right)}}} \right\}}} \right\rbrack} = 0} \\{{\sum\limits_{i = 1}^{N}\left\lbrack {{{- \frac{I_{i}}{R}}\left( {y - y^{i}} \right)} - {\frac{\left( {{\overset{\sim}{I}}_{i} - 1} \right)}{R}\left\{ {{\left( y^{i} \right)^{2}y} - {x^{i}{y^{i}\left( {x - x^{i}} \right)}}} \right\}}} \right\rbrack} = 0}\end{matrix};{\left( {x,y} \right) \in \mathcal{D}}} \right. & (156)\end{matrix}$

Equation 156 defines algorithm E which is a non linear, iterativealgorithm. The second approximation for the probability density functionf(x,y|z^(i)) given in section is f_(G) _(o) (x,y|z^(i)) in equation(149). This probability density function is obtained as a limit of thefirst Gaussian approximation, f_(G)(x,y|z^(i)), in the limitingcondition d_(o)→0. By using such expression, the problem simplifiesnoticeably; in fact the following partial derivatives result:$\begin{matrix}\left\{ \begin{matrix}{{\frac{\partial}{\partial x}\ln\quad{f\left( {x,\left. y \middle| z^{i} \right.} \right)}} = {- \frac{x - \mu_{x,i}}{I_{i0}}}} \\{{\frac{\partial}{\partial x}\ln\quad{f\left( {x,\left. y \middle| z^{i} \right.} \right)}} = {- \frac{y - \mu_{y,i}}{I_{i0}}}}\end{matrix} \right. & (157)\end{matrix}$where μ_(x,i)=x^(i), μ_(y,i)=y^(i) and I_(i0) defined in (148), dependson the i-th attenuation z^(i). With the result above, the MS locationestimate can be calculated in embodiments of the invention in closedform as follows: $\begin{matrix}{{\hat{x} = \frac{\sum\limits_{i = 1}^{N}\frac{x^{i}}{\mathcal{I}_{i\quad 0}}}{\sum\limits_{i = 1}^{N}\frac{1}{\mathcal{I}_{i\quad 0}}}};\quad{\hat{y} = \frac{\sum\limits_{i = 1}^{N}\frac{y^{i}}{\mathcal{I}_{i\quad 0}}}{\sum\limits_{i = 1}^{N}\frac{1}{\mathcal{I}_{i\quad 0}}}};\quad{\left( {\hat{x},\hat{y}} \right) \in \mathcal{D}}} & (158)\end{matrix}$

Equation 158 defines algorithm F which is a linear algorithm with aclosed form solution.

Extensions of the closed-form method

In the algorithm defined by equation (100) and used in embodiments ofthe present invention the estimated x (and respectively, y) MScoordinate is obtained as a weighted average of the x (and respectively,y) BTS coordinates, the signal of which is received by the MS. Thisalgorithm can be extended by means of the following formulation:$\begin{matrix}{{\hat{x} = \frac{\sum\limits_{i = 1}^{N}{w^{i}x^{i}}}{\sum\limits_{i = 1}^{N}w^{i}}};\quad{\hat{y} = \frac{\sum\limits_{i = 1}^{N}{w^{i}y^{i}}}{\sum\limits_{i = 1}^{N}w^{i}}};\quad{\left( {\hat{x},\hat{y}} \right) \in \mathcal{D}}} & (159)\end{matrix}$where w¹, . . . ,w^(N) are suitable weights assigned to each one of theN BTS's involved.

The weights in the algorithm (158) are calculated as the reciprocals ofthe terms I₁₀, . . . ,I_(N0). The weights being calculated in such a waythat the weights decrease as the attenuation increases, as can be seenalso from FIG. 11. It is possible to suggest that the signalstransmitted by BTS's closer to the MS undergo lower attenuation thansignals transmitted by BTS's located further away from the mobiledevice. This suggestion can be accounted for by imposing the rule thatthe weights in the generalized closed-form algorithm (159) are definedin such a way that, if the attenuation of the signal received from thei-th BTS is low, then w^(i) is high, and vice versa: $\begin{matrix}\left. \left. z^{i}\nearrow \right.\Rightarrow\left. w^{i}\searrow \right. \right. & (160)\end{matrix}$

Hereafter follows three empirical definitions for the weights w¹, . . .,w^(N) according to criterion suggested above as used in furtherembodiments of the present inventions.

1. Definition of Weights

The first definition of the weights which reflects the above rulescomprises defining the i-th weight as the reciprocal of the attenuationexperienced by the signal transmitted by the i-th BTS:w ^(i)=1/z ^(i)   (161)

This definition of weights is further enhanced in further embodiments ofthe present invention by the introduction of auxiliary variableparameters to be determined by using experimental measurements.

2. Definition of Weights

A further definition for the i-th weight used in embodiments of thepresent invention is to use the inverse of the estimated distance,{circumflex over (d)}^(i), between the MS and the i-th BTS, obtainedfrom the i-th observed attenuation z^(i):w ^(i)=1/d ^(i)   (162)

The level of the signal received by a MS {and the attenuation as well)does not only depend on the distance between MS and BTS, but also on thegain of the transmit antenna in the direction if the MS. However, if theN BTS's involved in the location calculation are omni-directional, thecontribution of the BTS antennas can be neglected and the attenuationcan be considered to be a function of the MS-BTS distance only. For thisapproximation, the i-th attenuation Z^(i) is equal to the i-th path-lossPL^(i)(d^(i)), and is dependent only on the distance between MS and thei-th BTS d^(i)={square root}{square root over ((x−x^(i))²+(y−y^(i))²)}:z ^(i) =PL ^(i)(d ^(i))   (163)

Neglecting the slow-fading in model (128), the following expression forthe weights as given in (162) produces the following: $\begin{matrix}{w^{i} = {\frac{1}{{\hat{d}}^{i}} = {10^{- \frac{{PL}^{i} - A}{B}} = 10^{- \frac{P_{t,\max}^{i} - P_{r}^{i} - A}{B}}}}} & (164)\end{matrix}$

A further alternative can be found when the weights defined in (164) areadjusted to fit experimental measurements. Two observations, describedbelow, can be made in order to produce more accurate weights to befound:

Since the MS coordinates are estimated using (159) with weights given by(162), the absolute values of the distances are not strictly needed. Inthe algorithm defined in equation (159) only the ratios 1/{circumflexover (d)}^(i)/Σ_(i=1) ^(N)/1/{circumflex over (d)}^(i) contribute to thelocation estimate. This means that if the same relative error is made inestimating each of the distances used in the location calculation, theresulting location estimate is not affected. This automatic errorcancellation factor means that going from outdoor to indoor or from openarea in a city to narrow city canyon should not greatly affect thelocation accuracy. This error cancellation factor improves the accuracyof the location estimate and minimizes the need to use a highlyoptimised and parameterised path-loss model in embodiments of theinvention.

For a location estimate calculation using the algorithm (159) andweights (164) only the functional behaviour of the distance as afunction of the received signal strength is needed. In a simplifiedexample, if the base station antenna height, central frequency, maximumradiated power are the same for each BTS used in the locationcalculation almost all the terms of the Okumura-Hata path loss modelcancel out and the resulting weights to be used in the locationcalculation simplify so that the weights can be found by the equation(165):w ^(i)=10^(−−P) ^(r) ^(i) ^(/B)   (165)3. Definition of Weights

A further embodiment of the present invention determines experimentallya functional relation between the weights, w^(i), and attenuation,z^(i). It is possible to collect a sufficient amount of received level(RXLEV) measurements, with corresponding GPS coordinates to be used toreference the exact MS locations, and from this information solve thereverse location problem: if the MS coordinates are known in equation(159), the weights w¹, . . . ,w^(N) can be determined as a function ofthe attenuation by minimizing the error between estimated MS coordinatesand exact MS coordinates. This approach fails where experimentalmeasurements are not available.

The location methods described above apply to a network made fromomni-directional cells, where the BTS's transmit antennas radiateisotropically in all directions. In sector cells, the assumptions madein an all omni-directional cell environment hold with a reasonable levelof approximation error only when the MS is located in the main loberegion (MLR). In other words when the mobile device is in the regionilluminated by the main beam of the BTS (directional) antenna. A sectorcell equipped with a transmit antenna having Half Power Beam Width(HPBW) of 60° has a MLR which approximately extends angularly ±60degrees about the direction towards which the BTS antenna is oriented.

In further embodiments of the present invention, the algorithms definedin equations (154), (156), (158), and (159) are applied if the N BTS'sinvolved in the location calculation, out of the total number of valuesreceived by the MS, are those that radiate the MS with the main beam oftheir transmit antennas.

Unfortunately it is not possible to determine whether the MS is in theMLR of a certain BTS from a single BTS received value because, in orderto obtain this information, the location of the MS is needed. However,in embodiments of the present invention it is possible to determine witha low error rate if the MS is in the MLR of a certain sector cell whenthe MS measures and compares the received levels (RXLEVs) from differentco-located sector cells. In this case, as the BTS antennas of co-locatedsector cells have different orientations, the MS is most likely locatedin the MLR of the cell the signal of which experiences the lowestattenuation among the co-located sector cells.

With the knowledge gathered by the method above the approximation foundfrom the use of omni-directional cells is accurate (e.g., andembodiments of the present invention may apply the algorithms defined inequations (154), (156), (158), and (159)) if the N BTS's involved in thelocation calculation are selected out of the ones measured by the MSaccording to the following procedure:

-   -   1. Select all omni-directional cells measured by the MS.    -   2. Select all sector cells, where no co-located cells are        measured by the MS at the same time.    -   3. Select the sector cell with the lowest attenuation out of the        co-located sector cells, where co-located cells are detected by        the MS.

If one or more of the discarded co-located cells produce a receivedlevel value attenuation that differs from the selected cell by an amountby a certain range (say 5 dB), select these co-located cells also andthe average of the attenuation in embodiments of the present inventionor in further embodiments of the present invention use the attenuationvalues of all of the selected co-located cells.

CI and CI+TA Based Algorithms

These algorithms will be described in more detail hereinafter. Thesealgorithms are linear and rely only on serving cell information. Theyboth deliver a location estimate and a confidence region for thelocation estimate.

These algorithms are considered to act upon a mobile device beingserviced by one generalised sector cell with a coverage or hearabilityarea defined analytically or by practical measurements. The generalisedcell is a simplification of the effect produced by the BTS antennaangular gain to both transmit and receive as shown by the practicalantenna gain shown in FIG. 5.

As explained previously and as shown in FIG. 10 the generalised cell isdefined by a series of parameters comprising, φ the sector orientationin degrees from the x-axis, Δφ the sector angular width, R_(F) thesector front radius, and R_(B) the sector back radius.

The parameters defining the generalised sector may be found fromcoverage maps or coverage prediction tools. The input informationrequired by the method comprises, the Cell Identity of the serving BTS(CI), the serving BTS coordinates; and the map of the serving cell.

The steps performed to calculate the generalised cell parameters from acoverage map of the cell comprise:

-   -   1. Selecting the Coordinates of the Serving BTS Using the CI:        x_(s), y_(s).

Given the CI of the serving cell, the x-y coordinates of the cell areidentified the CI, (x_(s),y_(s)), and are retrieved from the database ofBTS-coordinates.

2. Selecting the Serving Cell Map: S

Analogously as in the previous step, the CI is used to select theserving map, S, corresponding to the current CI. S is typicallydetermined by a network planning tool which uses a large set ofinformation such as BTS configuration parameters, 3-dimensional terrainmaps, tuned propagation models, etc. and, in order to determine theserving area, takes into account also the presence of other BTS's in thenetwork.

The serving map represents the geographical region where the cellidentified by CI is serving. In an arbitrary x-y Cartesian system, S canbe obtained by dividing the region into elements of area (Δx)×(Δy) andrepresenting each element by the coordinates of its centre:S:{x_(n),y_(n)}; n=1, . . . ,N_(S)   (166)

The coordinates (x_(n),y_(n)) represent the centre of the n-th pixel ofarea (Δx)×(Δy) where the cell is serving.

An example of serving area is shown in FIG. 13. FIG. 13 shows a basetransceiver station (BTS) 1301, a main coverage area 1303 and minorcoverage areas 1305. The base transceiver station 1301 lies within themain coverage area 1303. Smaller coverage areas 1305 are positionedadjacent to but not touching the edges of the main coverage area 1303.All of the coverage areas are divided into a plurality of elements 1309by a grid system 1307. The total of all of the coverage areas dependenton the BTS is known as the serving area for the BTS.

3. Determination of the Coordinates of the Centre of Mass of the ServingMap: x_(mc),y_(mc)

The (x,y) coordinates of the serving cell's mass centre are formallydefined as follows: $\begin{matrix}{{x_{MC} = {\frac{1}{\mathcal{M}(\mathcal{S})}{\int_{\mathcal{S}}{x{\mathbb{d}x}{\mathbb{d}y}}}}};\quad{y_{MC} = {\frac{1}{\mathcal{M}(\mathcal{S})}{\int_{\mathcal{S}}{y{\mathbb{d}x}{\mathbb{d}y}}}}}} & (167)\end{matrix}$where M(S) is the area of the serving map:M(S)=∫_(S) dxdy   (168)

By using the divided elements defined in equation (166) the coordinatesof the mass centre can be calculated as follows: $\begin{matrix}{{x_{MC} \simeq {\frac{1}{N_{S}}{\sum\limits_{n = 1}^{N_{S}}x_{n}}}};\quad{y_{MC} \simeq {\frac{1}{N_{S}}{\sum\limits_{n = 1}^{N_{S}}y_{n}}}}} & (169)\end{matrix}$

4. It is much easier to calculate the desired estimates in a polarreference system (ρ_(i),θ_(i)) originated in the serving BTS ofcoordinates (x_(s),y_(s)). In this polar coordinate system ρ_(i) is thedistance of the point (x_(i),y_(i)) from the BTS and θ_(i) is the anglemeasured counter-clockwise from the x-axis: $\begin{matrix}{{\rho_{i} = \sqrt{\left( {x_{S} - x_{i}} \right)^{2} + \left( {y_{S} - y_{i}} \right)^{2}}};{\theta_{i} = {\tan^{{- 1}\frac{y_{S} - y_{i}}{x_{S} - x_{i}}}.}}} & (170)\end{matrix}$

In the polar reference system the serving map is represented by a set ofpoints:S_(polar):{ρ_(n),θ_(n)}; n=1, . . . ,N_(S).   (71)

5. Determination of the Main Direction for the Entire Serving Map:φ_(mc)

The main direction θ_(mc) provides an indication of the serving sector'sbearing, φ_(s). This angle is used as a reference direction fornormalizing the angular coordinates of the cell's coverage areaelements. θ_(mc) can be approximated as the orientation of the cell'smass centre from the BTS coordinates.

In the polar reference system (ρ, θ) the mass centre of the serving areahas coordinates (ρ_(mc),θ_(mc)) such that: $\begin{matrix}{{\rho_{MC} = \sqrt{\left( {x_{S} - x_{MC}} \right)^{2} + \left( {y_{S} - y_{MC}} \right)^{2}}};{\theta_{MC} = \tan^{{- 1}\frac{y_{S} - y_{MC}}{x_{S} - x_{MC}}}}} & (172)\end{matrix}$

6. To further simplify the actual calculation of the estimates the newpolar reference system is rotated in such way that zero and 2π anglesare as far as possible from the most important serving area. After thisrotation all the needed estimates can be calculated using a known simplesorting algorithm.

The rotation of the m-th pixel with angular coordinate θ_(m), is definedas follows: $\begin{matrix}{\delta_{m} = \left\{ {\begin{matrix}{\theta_{m} - \theta_{R} + {2\pi}} & {{{{if}\quad\theta_{m}} - \theta_{R}} < 0} \\{\theta_{m} - \theta_{R}} & {{{{if}\quad\theta_{m}} - \theta_{R}} \geq 0}\end{matrix}{where}} \right.} & (173) \\{\delta_{R} = \left\{ \begin{matrix}{\theta_{MC} + \pi} & {{{if}\quad\theta_{MC}} < \pi} \\{\theta_{MC} - \pi} & {{{if}\quad\theta_{MC}} \geq \pi}\end{matrix} \right.} & (174)\end{matrix}$

The resulting δ_(m)'s are in the range 0<δ_(m)<2π. The directionsδ_(m)'s closest to the main direction θ_(mc) are associated to δ_(m)≅πand the directions furthest away from θ_(mc) are associated to δ_(m)≅0and δ_(m)≅2π.

Next step is to sort the points in order of increasing distance from theserving BTS:{ρ_(l),δ_(l)}_(t=1) ^(N) ^(s) =sort {ρ_(l)}.   (175)

At this position the serving map is represented by an ordered set ofpoints:S_(polar): {ρ_(n),δ_(n)}; n=1, . . . ,N_(S).   (176)7. Determination of the Front Radius, R_(F):R _(F)=ρ _(l|l=┌NS┐) ,   (177)where ┌•┐ rounds • to the nearest integer towards plus infinity fusingthe known ceiling operation where all numbers are rounded up towards thenext positive integer) and, γ′ determines the desired fraction of pointsin set S_(polar) to have the distance from serving BTS less than R_(F).A good value for the parameter, γ′ is between 0.95 and 0.98, in someembodiments of the invention.

8. Determination of back radius (R_(B)), orientation as a function ofdistance from serving BTS (φ_(s)(d)), and angular width as a function ofdistance from serving BTS (Δφ_(s)(d)).

First it is helpful to define the parameters of the circular crowndepicted in FIG. 14. FIG. 14 shows a constant timing advance (TA) regiondefined by the difference in area between two concentric circles with acommon origin 1401. The first circle 1403 having a radius R₁ 1407 andthe second circle 1405 having a radius R₁+R₂.

The mathematical definition for the crown in the polar reference systemis:C: {(ρ,θ)εR ² :R _(inf) ≦ρ≦R _(sup)}  (178)Where R_(inf) is inner radius and R_(sup) is outer radius. In theexample shown in FIG. 14 the first circle radius 1407 R₁=R_(inf), andthe second circle radius R₁+R₂=R_(sup).

This crown is therefore represented by a set of points:S _(crown): {ρ_(n),δ_(n) }; n=1, . . . ,N _(crown)(N _(crown) ≦N _(S)).  (179)

9. To determine the back radius (R_(B)) the following steps arefollowed:

-   -   (a) set i=0    -   (b) The number of points in the crown defined in equation (179)        is calculated using R_(inf)=i·Δ and R_(sup)=(i+2)·Δ, where Δ is        the distance of digitized elements in x and y directions.    -   (c) If the serving area is omni-directional at the distance        R_(sup), there are approximately N_(ideal)=π(R_(sup) ²−R_(inf)        ²)/Δ² points in the crown. If the actual number of points in        S_(crown), N_(crown) is less than γ′N_(ideal,), the serving area        is no longer considered to be omni-directional. The value of        γ′<1.0 and is based on experience a good value being 0.75. If        the serving area is not considered to be omni-directional the        back radius is chosen to be R_(B)=R_(sup).    -   (d) If the serving area can be considered omni-directional at        distance R_(sup), set i=i+1 and go to (b), otherwise the        estimate of Back radius is ready.

10. To determine the orientation, φ_(s)(d), the angles in the circularcrown are sorted in order of either increasing or decreasing angle;{ρ_(l),δ_(l)}_(i=1) ^(N) ^(crown) =sort {δ_(l)},   (180)

From the definition of the crown R_(inf)=d−2 Δ and R_(sup)=d+2 Δ. Theestimate for the φ_(s)(d) is therefore the median value of the anglesand it is obtained as:φ_(S)(d)=δ_(l|l=┌(N) _(crown) _(/2┐).   (181)

11. The Angular width as a function of distance from serving BTS,Δφ_(s)(d) is calculated using the same sorted set of points in thecircular crown as is used for φ_(s)(d). Angular width, Δφ_(s)(d), iscalculated as follows:Δφ_(S)(d)=max{∥δ₁−φ_(S)(d)∥,∥δ_(N) _(crown) −φ_(S)(d)∥}.   (182)

12. As the last step it is necessary to rotate φ_(s)(d) and Δφ_(s)(d)back to the original co-ordinate reference system. This is done bysimply adding θ_(R) to the obtained values.

Where coverage maps are not available, it is possible to analyticallyestimate the cell front radius, R_(F), and cell back radius, R_(B), withthe method presented below.

As described earlier, such as by equation (33), the average receivedpower, P_(R), by a mobile device a distance, d, from the serving BTS canbe expressed as:P _(R)(d)=P _(T) +G′−PL(d)   (183)

The BTS transmits a known transmitted power P_(T) and the gain of theantenna installed at the serving BTS in the direction of the MS is G.PL(d) is the path-loss affecting the signal power as it propagates fromserving BTS to MS. The accepted model, such as those used in equations(126) to (129) above, assumes the path-loss increases logarithmicallywith the distance d:PL(d)=A+B log d   (184)

The sector front radius can be defined as the distance at which theaverage power received by a MS, P_(R), is above a certain thresholddefining the cell edges, P_(R) ^(th). The effect of potential shadowfading is allowed for by including a shadow fade margin FM_(σ) (allquantities are in logarithmic units)P _(R) =P _(R) ^(th) +FM _(σ)  (185)

Under the assumption of log-normal slow fading, the fade margin can bedefined as FM_(σ)=zσ, where σ is the standard deviation of the slowfading and z is such that${F(z)} = {{1 - {Q(z)}} = {1 - {\frac{1}{\sqrt{2\pi}}{\int_{z}^{\infty}{{\mathbb{e}}^{{- x^{2}}/2}{\mathbb{d}x}}}}}}$is the radius estimate's reliability.

According to equations (183) and (185) the front cell radius inkilometres can be then expressed as $\begin{matrix}{R_{F} = 10^{\frac{F_{T} - F_{R}^{th} - {F\quad M_{c}} + G_{m} - A}{B}}} & (186)\end{matrix}$where G_(m) is the maximum antenna gain in dB.

The sector back radius can be calculated from the front radius byincluding information on the radiation pattern of the antenna installedat the serving BTS site. FIG. 5 shows an annotated version of FIG. 4, aplot of antenna gain in decibels against angle measured in degrees. Theantenna gain, symmetrical about the angle θ=0, comprises a main lobe501, and four side lobes 503, 505, 507, 509. The main lobe 501 has amaximum value of G_(m) at θ=0, whilst none of the side lobes 503, 505,507, 509 have gains greater than ρG_(m). ρ is the maximum Back-to-FrontRatio, that is the maximum ratio between the antenna gain in thedirection of maximum radiation, G_(m) and the gain of the antenna in thedirections outside of the main lobe. It is evident that the cell frontradius, R_(F), is related to the maximum gain in the main lobedirections, G_(m), and the cell back radius, R_(B), is related to themaximum gain in the directions outside of the main lobe ρG_(m).

By using the simple propagation model (183) and the following roughapproximation of the antenna gain G(θ) in dB (G_(m)=10 log(g_(m)))$\begin{matrix}{{G(\theta)} \simeq \left\{ \begin{matrix}{G_{m};} & {0 \leq \theta < \theta_{3{dB}}} \\{{G_{m} + {10{\log(\rho)}}};} & {\theta_{3{dB}} \leq \theta \leq {{2\pi} - \theta_{3{dB}}}} \\{G_{m};} & {{{2\pi} - \theta_{3{dB}}} < \theta \leq {2\pi}}\end{matrix} \right.} & (187)\end{matrix}$it is possible to write the following upper limit for R_(B):R _(B) ≦ρ×R _(F)   (188)

Typical values used in GSM900 are P_(T)=50 dBm, P_(R) ^(th)=−95 dBm. Fora BTS height of 30 meters, typical values for variables are A=124.5 andB=35.7. In the mobile radio environment, a is often assumed equal to 8dB. For this example the BTS antenna has HPBW 65 degrees, Back-to-FrontRatio −18 dB (ρ=0.0158) and maximum gain 12 dB. With these values andfor z=0.675 (i.e., F(z)=0.75 equivalent to 75% cell radius reliability)formulas (186) and (188) give R_(F)=5.7425 kilometers and R_(B)≅90meters.

From the examples shown, it is evident that the representation of thesector as shown in FIG. 10 might approximate poorly the real cellcoverage. In particular, Δφ_(s) is usually too large to represent thesector width in those regions far away from the BTS coordinates. Toimprove this inaccuracy in estimation, it is possible to define theserving sector width as a function of the distance from the serving BTS,Δφ_(s)(d). By taking this into account it is possible to neglect in thelocation algorithm those regions which clearly lie outside of theserving area. The same applies for the sector direction, φ_(s), whichcan be assumed to be a function of the distance from the BTS site,φ_(s)(d).

R_(F), R_(B), φ_(s) and φ_(s)(d) (eventually dependent on the distance,d) create a analytical estimation of the simplified borders of theserving cell as follows: $\begin{matrix}{{\mathcal{S}\left( {x,y} \right)}\text{:}\quad\left\{ \begin{matrix}{{{d\left( {x,y} \right)} = R_{F}};} & {0 \leq {{{\psi\left( {x,y} \right)} - \phi_{S}}} \leq {\Delta\phi}_{S}} \\{{{d\left( {x,y} \right)} = R_{B}};} & {{{{\psi\left( {x,y} \right)} - \phi_{S}}} > {\Delta\phi}_{S}}\end{matrix} \right.} & (189)\end{matrix}$

In embodiments of the present invention the above estimations of thegeneralised cell parameters are used to estimate the location of themobile device (MS). The algorithm used in embodiments of the presentinvention is referred to as a Classic CI-TA Location Algorithm, asopposed to a Map-Aided CI-TA Location Algorithm.

As explained previously the input parameters to the algorithm used inembodiments of the present invention are CI, TA (if available) andconfidence coefficient of the location estimate, ξ. As also explainedpreviously the output of the calculation is a location estimate and aconfidence region (e.g., a geographical region within which the real MSlocation is situated within a degree of confidence ξ). The result of thelocation calculation differs if the cell is sectorized oromni-directional.

Embodiments of the present invention described below are capable ofproviding these results when both CI and TA are available, furtherembodiments also described below provide results when only CI isavailable.

If the serving cell is sectorized, the location estimate is calculatedby combining the serving BTS coordinates, x_(s) and y_(s), with anestimate of the MS-to-serving BTS distance, {circumflex over (d)}, andan estimate of the angular coordinate of the MS from the serving BTSsite, {circumflex over (ψ)}: $\begin{matrix}\left\{ \begin{matrix}{\hat{x} = {x_{S} + {\hat{d}\quad\cos\quad\hat{\psi}}}} \\{\hat{y} = {y_{S} + {\hat{d}\quad\sin\quad\hat{\psi}}}}\end{matrix} \right. & (190)\end{matrix}${circumflex over (d)} and {circumflex over (ψ)} are the estimateddistance and angle respectively as described below.

Given the TA information, the method described previously can be used todetermine an estimate of the distance between MS and serving BTS and theradii of a circular crown (C) centred at the BTS coordinates where themobile station can be located with a confidence γ.

The distance estimate, {circumflex over (d)}, is calculated as the 50-thpercentile (or median value) of the real distance, d, i.e. {circumflexover (d)} is such thatPr(d≦{circumflex over (d)})=½  (191)

Other analogous definitions for {circumflex over (d)}, such as the meanof d, are possible. Using the results described in European PatentApplication number 102251 which is hereby incorporated by reference, themedian estimated distance is:{circumflex over (d)}=d _(TA) +TA _(c) ; TA _(c) =−X _(1/2)   (192)where X_(1/2) is the median value of the TA measurement error,X=d_(TA)−d, and $\begin{matrix}{d_{TA} = \left\{ {\begin{matrix}{{\frac{1}{4} \times \left( {{cT}_{b}/2} \right)} \simeq {138m}} & {{{if}\quad{TA}} = 0} \\{{TA} \times \left( {{cT}_{b}/2} \right)} & {{{if}\quad{TA}} > 0}\end{matrix},} \right.} & (193)\end{matrix}$

-   -   T_(b)=3.69 μs is the bit period and c=3×10⁸ m/s is the speed of        light.

In absence of any other information, the angular coordinate of the MScan be estimated with the orientation of the sector, i.e. {circumflexover (ψ)} can be set equal to φ_(s).

If the sector orientation is a function of the distance from the servingBTS, then {circumflex over (ψ)} is the sector orientation at a distanceequal to the one estimated on the basis of the TA, i.e. {circumflex over(ψ)}=φ_(s)({circumflex over (d)}).

Further means for estimating the MS angular coordinate can beincorporated into the algorithm found in embodiments of the presentinvention; for example, {circumflex over (ψ)} can in further embodimentsof the invention be determined by processing signal level measurements(RXLEVs) performed by the MS or, in the future, by using theangle-of-arrival information made available by (smart) antenna arraysinstalled at the serving BTS site.

Associated to the estimated MS location is the confidence region R asshown in FIG. 8 and comprising the following parameters: origin locatedat point with coordinates (x₀,y₀); inner radius, R₁; uncertainty radius,R₂; orientation angle measured counter-clockwise from x axis, α, andinclusion angle defining the width of the sector, β.

Analogously to the distance estimate, the confidence region can bedetermined by using the method described in European Patent Applicationnumber 102251 which is hereby incorporated by reference.

The origin of the confidence region is at the serving BTS site:$\begin{matrix}\left\{ \begin{matrix}{x_{0} = x_{S}} \\{y_{0} = y_{S}}\end{matrix} \right. & (194)\end{matrix}$

Given the measured TA and the confidence coefficient for the distanceestimate γ, the confidence interval of the distance estimate({circumflex over (d)}) can be determined from the statisticalproperties of the TA measurement provided by the known map aided CI-TAestimation technique. The confidence interval for the distance estimateis determined by r_(i) and r_(s) defined as $\begin{matrix}\left\{ \begin{matrix}{r_{i} = {{TA}_{c} + X_{{({1 + \gamma})}/2}}} \\{r_{s} = {{- {TA}_{c}} - X_{{({1 - \gamma})}/2}}}\end{matrix} \right. & (195)\end{matrix}$where X_((1±γ)/2) is the 100 1±γ/2% percentile of the TA measurementerror. r_(i) and r_(s) are such that the true MS-serving BTS distance,d, falls within the confidence interval [{circumflex over (d)}−r_(i),{circumflex over (d)}+r_(s)] with a probability, γ as shown in FIG. 15):Pr({circumflex over (d)}−r _(i) ≦d≦{circumflex over (d)}+r_(s))=γ  (196)

FIG. 15 shows a base transceiver station (BTS) 1501, a location estimate{circumflex over (d)} 1507 which lies within the area defined by twoarcs 1503, 1505. The larger arc has a radius d_(sup), the smaller archas a radius d_(inf).

In embodiments of the present invention the confidence region parameters(R₁, R₂, α, and β) are calculated slightly differently depending if theserving cell back radius (R_(B)) is smaller or larger than the lowerlimit of the confidence interval for the distance estimate,d_(INF)={circumflex over (d)}−r_(i).

The cell back radius R_(B) is not strictly needed, but provides enoughinformation to improve the reliability of the method. Therefore infurther embodiments of the present invention the value of the backradius R_(B) is set to zero.

If d_(INF) is larger than the serving cell back radius, R_(B), the MS issufficiently far away from the serving BTS and is probably locatedoutside of the back radius region. In this case, only the portion ofcell in the direction of main radiation lobe of the serving BTS antennais included in the confidence region.

The inner radius of the confidence region is set to R₁={circumflex over(d)}−r_(i). However R₁={circumflex over (d)}−r_(i)=d_(TA)−X_((1±γ)/2)can be negative (i.e., if γ≈1 and the TA error statistics have longtails X_((1±γ)/2) can be larger than d_(TA)), and if R₁ is calculated tobe negative, R₁ should be set to zero. Moreover, when TA=0 it issensible to define the confidence region as a circle centred on theserving BTS coordinates; in which case R₁ is set to zero. In summary, R₁is defined as follows: $\begin{matrix}{R_{1} = \left\{ \begin{matrix}{0} & {{{if}\quad{TA}} = 0} \\{\max\left\{ {0,{\hat{d} - r_{i}}} \right\rbrack} & {{{if}\quad{TA}} > 0}\end{matrix} \right.} & (197)\end{matrix}$

The uncertainty radius of the confidence region (R₂) is set equal to thewidth of the confidence interval of the distance estimate, r_(i)+r_(s).However, as a consequence of the adjustments made in (197) to determineR₁, the correct definition for R₂ isR ₂ ={circumflex over (d)}+r _(s) −R ₁.   (198)

The orientation angle of the confidence region (α) is determined by theorientation of the serving sector and the sector angular width. Thewidth of the confidence region β is twice the serving sector width:$\begin{matrix}\left\{ {\begin{matrix}{\alpha = {{\phi_{S}\left( \hat{d} \right)} - \overset{\_}{{\Delta\phi}_{S}}}} \\{\beta = {2\overset{\_}{{\Delta\phi}_{S}}}}\end{matrix},} \right. & (199)\end{matrix}$where {overscore (Δφ_(s))} is defined to take into accountconservatively the eventual dependence of the cell width on the distancefrom the serving BTS, d: $\begin{matrix}{\overset{\_}{\Delta\quad\phi_{S}} = {\max\limits_{R_{1} \leq d \leq {R_{1} + R_{2}}}{\Delta\quad{\phi_{s}(d)}}}} & (200)\end{matrix}$

The confidence region has an angular width equal to twice the angularwidth of the serving cell. The sectorization of the cell is thereforetaken into account. The underlying reason for this is the assumptionthat all the mobile devices in communication with the BTS of interestare located in an arc defined angularly −Δφ_(s)≦ψ(x,y)−φ_(s)({circumflexover (d)})≦Δφ_(s) (with an exception covered by equation (189) whenR_(B)=0). This also produces the result that the confidence coefficientused to determine the circular crown, γ, is equal to the confidencecoefficient ξ used to define distances R₁ and R₂.

If d_(INF) is smaller than the serving cell back radius, R_(B), there isa large probability that the actual location of the MS lies in the backradius region. In this case the location estimate is calculated aspreviously but the confidence region is defined in such a way to includethe whole back radius region. This means that the width of theconfidence region β is 2π and the orientation angle α provides redundantinformation as the confidence arc is defined as the whole circle (α isset to zero). Furthermore, the inner and the uncertainty radii are givenas follows: $\begin{matrix}\left\{ \begin{matrix}{R_{1} = 0} \\{R_{2} = {\hat{d} + {r_{s}.}}}\end{matrix} \right. & (201)\end{matrix}$

In practical terms, the definitions above simply state that when the MSis very close to the serving BTS site, the effects of the back radiusregion are considered by treating the serving cell as if it wasomni-directional.

If the serving cell is omni-directional the concepts of sectororientation and angular width are meaningless. The best locationestimate is provided by using the co-ordinates given by the serving BTSsite: $\begin{matrix}\left\{ \begin{matrix}{\hat{x} = x_{S}} \\{\hat{y} = y_{S}}\end{matrix} \right. & (202)\end{matrix}$where the confidence region is defined by a circle with radius equal tothe upper limit of the confidence interval of the distance estimate:$\begin{matrix}\left\{ {\begin{matrix}{x_{0} = x_{S}} \\{y_{0} = y_{S}}\end{matrix};\left\{ {\begin{matrix}{\alpha = {{Any}\quad{value}}} \\{\beta = {2\pi}}\end{matrix};\left\{ \begin{matrix}{R_{1} = 0} \\{R_{2} = {\hat{d} + r_{s}}}\end{matrix} \right.} \right.} \right. & (203)\end{matrix}${circumflex over (d)} and r_(s) are calculated as in the case ofsectorized cell using equations (192) and (195).

The calculations described above make use of statistical informationfrom the TA measurement error. If the TA measurement error statisticalinformation is not available, the definitions for location estimate andconfidence region given above still apply; however the confidencecoefficient is meaningless and {circumflex over (d)}, r_(i) and r_(s)can be defined only on the basis of the TA quantization rule:$\begin{matrix}{{\hat{d} = d_{TA}}{r_{i} = {r_{s} = \left\{ \begin{matrix}{{{\frac{1}{4} \times \left( {{cT}_{b}/2} \right)} \cong {138m\quad{if}\quad{TA}}} = 0} \\{{\frac{1}{2} \times \left( {{cT}_{b}/2} \right)} \cong {277m\quad{if}\quad{TA}} > 0}\end{matrix} \right.}}} & (204)\end{matrix}$

The algorithm presented above estimates the MS coordinates and theconfidence region when both CI and TA are available. If no TAinformation is available, the MS location estimate and its confidenceregion can still be determined by using the information carried by theCI only. In particular, the confidence region can be defined if anestimate on the cell radius, R_(F) is available while the locationestimate can be provided even if R_(F) is not known.

Where there is no TA information available, embodiments of the presentinvention examine the serving BTS to determine whether the serving cellis sectorized or not:

If the cell is omnidirectional and the TA is not available the locationestimate is according to embodiments of the invention taken to be at theBTS coordinates: $\begin{matrix}\left\{ \begin{matrix}{\hat{x} = x_{S}} \\{\hat{y} = y_{S}}\end{matrix} \right. & (205)\end{matrix}$

The confidence region of this estimate is a circle centred on the BTSsite. The radius of such region is obtained by scaling R_(F) by a factor{square root}{square root over (ξ)}. The factor is chosen in order thata fraction ξ of the total area of the circular cell of radius R_(F) isincluded in the confidence region: $\begin{matrix}\left\{ \begin{matrix}{{x_{0} = x_{S}};} & {y_{0} = y_{S}} \\{{\alpha = {{Any}\quad{value}}};} & {\beta = {2\pi}} \\{{R_{1} = 0};} & {R_{2} = {\sqrt{\xi}R_{F}}}\end{matrix} \right. & (206)\end{matrix}$

If the cell is sectorized and the TA is not available, two alternativesexist for the location estimate.

The first alternate as used in embodiments of the invention select theMS location estimate to be at the serving BTS site. In this estimate nocell radius (R_(F)) information is needed. $\begin{matrix}\left\{ \begin{matrix}{\hat{x} = x_{S}} \\{\hat{y} = y_{S}}\end{matrix} \right. & (207)\end{matrix}$

If the R_(F) is given, embodiments of the present invention calculatethe estimated MS location as the mass centre of a simplified cell withthe same shape as represented in FIG. 10 but with the back radius andthe front radius scaled by a factor {square root}{square root over (ξ)},to assure that only a fraction ξ of the total area of the original cellis considered.

The coordinates of the serving cell's mass centre are defined asfollows: $\begin{matrix}{{x_{MC} = {\frac{1}{M(S)}{\int_{S}{x\quad{\mathbb{d}x}{\mathbb{d}y}}}}};{y_{MC} = {\frac{1}{M(S)}{\int_{S}{y\quad{\mathbb{d}x}{\mathbb{d}y}}}}}} & (208)\end{matrix}$where S is the border of the home cell, as explained in equation (189),and M(S) its area. By neglecting the dependence of the sectororientation and width on the distance and assuming a constant sector'sangular width {overscore (Δφ_(s))} defined as $\begin{matrix}{\overset{\_}{{\Delta\phi}_{s}} = {\max\limits_{R_{B} \leq d \leq R_{F}}{{\Delta\phi}_{s}(d)}}} & (209)\end{matrix}$

it can be shown that M(S)=R_(F) ²{overscore (Δφ_(s))}−(π−{overscore(Δφ_(s))})R_(B) ² and the integrals in (208) solved in polar coordinatesgive the following expressions for the location estimate {circumflexover (x)}=x_(mc) and {circumflex over (γ)}=y_(mc) $\begin{matrix}\left\{ \begin{matrix}{\hat{x} = {x_{S} + {\frac{2}{3}\frac{\left( {R_{F}^{3} - R_{B}^{3}} \right)\sin\quad\overset{\_}{{\Delta\phi}_{S}}}{{R_{F}^{2}\overset{\_}{{\Delta\phi}_{S}}} - {\left( {\pi - \overset{\_}{{\Delta\phi}_{S}}} \right)R_{B}^{2}}}\cos\quad\phi_{S}}}} \\{\hat{y} = {y_{S} + {\frac{2}{3}\frac{\left( {R_{F}^{3} - R_{B}^{3}} \right)\sin\overset{\_}{\quad{\Delta\phi}_{S}}}{{R_{F}^{2}\overset{\_}{{\Delta\phi}_{S}}} - {\left( {\pi - \overset{\_}{{\Delta\phi}_{S}}} \right)R_{B}^{2}}}\sin\quad\phi_{S}}}}\end{matrix} \right. & (210)\end{matrix}$

The confidence region has the same shape as shown in FIG. 8, with theinner radius equal to zero (R₁=0). The origin of the confidence region(x_(o),y_(o)) is not at the BTS coordinates (x_(s),y_(s)) but is shiftedalong the axis defined by the sector orientation φ_(s) (in the directionopposite to the front side of the cell) by a distance R_(s)′ from theBTS coordinates. By using this definition of the confidence region, theback radius region is included (at least partially). The angular widthof the confidence region, β, is defined as twice the angleΔφ_(s)′≠Δφ_(s), which is calculated in light of the different locationof confidence region's origin and serving BTS coordinates. Theconfidence region's orientation, ox, is defined according to the cellorientation, φ_(s), and the new variable Δφ_(s)′: $\begin{matrix}\left\{ \begin{matrix}{{x_{0} = {x_{S} - {R_{B}^{\prime}\cos\quad\phi_{S}}}};} & {y_{0} = {y_{S} - {R_{B}^{\prime}\sin\quad\phi_{S}}}} \\{{\alpha = {\phi_{S} - {\Delta\phi}_{S}^{\prime}}};} & {\beta = {2{\Delta\phi}_{S}^{\prime}}} \\{{R_{1} = 0};} & {R_{2} = {\sqrt{\xi}\left( {R_{F} + R_{B}^{\prime}} \right)}}\end{matrix} \right. & (211)\end{matrix}$

-   -   R_(B)′ and Δφ_(s)′ are determined from the geometry as shown in        illustrated in FIG. 16.

FIG. 16 shows the geometry used in embodiments of the present inventionfor calculating the confidence region. The figure comprises a circle1603, a first triangle 1601, a second triangle 1605, and a circularsegment 1607. The circle has an origin, 1611, at the BTS location,x_(s),y_(s), and has a radius R_(B). The first comprises a first vertex1609 at location x₀,y₀, a first side 1623, connected at one end to thefirst vertex, 1609 which passes the origin of the circle, x_(s),y_(s), asecond side 1621, connected at one end to the first vertex and arrangedso that the angle at the first vertex 1609 defined by the two sides1623, 1621 is Δφ_(s)′. The vertex 1609 of the first triangle is locateda distance R_(B)′ from the origin of the circle. The second triangle1605 lies within the area of the first triangle 1601, and comprises afirst vertex located at the origin of the circle, a first side 1629connected at one end to the second triangle's first vertex and formingpart of the first triangle's first side, and a second side 1627connected at one end to the first vertex and arranged to have an anglebetween the first side 1629 and the second side 1627 of {overscore(Δφ_(s))}. Both the first triangle 1601 and the second triangle 1605,also comprise a common right angle vertex 1631, and a common side 1625.The common side being arranged as the line perpendicular to bothtriangles' first sides and connected to the point where the first andsecond triangle's sides intersect 1635. The segment 1607 comprises tworadii of length R_(f) with an origin at the origin of the circle 1603,x_(s),y_(s), and an arc defined between the ends of the radii. The firstradius is common to the second triangle's second side 1605, and thesecond radius is an angle 2{overscore (Δφ _(s) )}, from the first. Usingsimple trigonometry, the height of the first and second triangle'ssecond side can therefore be shown to be equal to R _(F) sin({overscore(Δφ_(s))}). FIG. 16 further has a line 1651 connected to andperpendicular to the first triangle's first side 1623 at the circle'sorigin 1611, and a second end at the intersection of the line and thefirst triangle's second side. The line 1651 has a length h.

Depending on whether {overscore (Δφ_(s))} is larger or smaller than π/2,embodiments of the present invention determine R_(B)′ and Δφ_(s)′ asexplained below.

If R_(F)>>R_(B) and {overscore (Δφ_(s))}>π/2, it can be seen that only avery small area of the circle defined by R_(B) is left out of theconfidence region if R_(B)′ is defined as being equal to R_(B):R_(B)′≡R_(B).   (212)

Using this definition for R_(B)′ the analytical expression for Δφ_(s)′is calculated from the geometry shown in FIG. 16. The result is:$\begin{matrix}{{\Delta\phi}_{S}^{\prime} = {\pi - {{\cos^{- 1}\left\lbrack \frac{{R_{F}{\cos\left( {\pi - \overset{\_}{{\Delta\phi}_{S}}} \right)}} - R_{B}^{\prime}}{\sqrt{\left( {{R_{F}{\cos\left( {\pi - \overset{\_}{{\Delta\phi}_{S}}} \right)}} - R_{B}^{\prime}} \right)^{2} - {R_{F}^{2}{\sin^{2}\left( {\pi - \overset{\_}{{\Delta\phi}_{S}}} \right)}}}} \right\rbrack}.}}} & (213)\end{matrix}$

Using the geometry shown in FIG. 16 Δφ_(s)′ is obtained as follows:$\begin{matrix}{{\Delta\phi}_{S}^{\prime} = {\tan^{- 1}\left\lbrack \frac{{R_{F}{\sin\left( \overset{\_}{{\Delta\phi}_{S}} \right)}} - h}{R_{F}{\cos\left( \overset{\_}{{\Delta\phi}_{S}} \right)}} \right\rbrack}} & (214)\end{matrix}$and using this result R_(B)′ is easy to calculate as follows:$\begin{matrix}{R_{B}^{\prime} = \frac{h}{\tan\left( \overset{\_}{{\Delta\phi}_{S}^{\prime}} \right)}} & (215)\end{matrix}$

The length h can be approximated as h≈R_(B). This approximation is onlyperformed after checking that the confidence region is too small anddoes not include a large enough part of the circle with radius of R_(B)(i.e. the value of R_(B)′ is too small) or conversely that R_(B)′ is toolarge leading to a too large confidence region. To check for these casestwo chosen parameters are introduced: δ_(min) (taking care of the toosmall R_(B)′ example—where 1.5 is a typical choice for δ_(min)) andδ_(max) (taking care of the too large R_(B)′ example—there 4-5 is atypical choice for δ_(max)).

If R_(B)′ turns out to be too small, i.e. if δ_(min)R_(B)>R_(B)′, thenR_(B)′ is re-defined in embodiments of the invention asR_(B)′≡δ_(min)R_(B)   (216)and Δφ_(s)′ is obtained using the geometry as shown in FIG. 16:$\begin{matrix}{{\Delta\phi}_{S}^{\prime} = {\tan^{- 1}\left\lbrack \frac{R_{F}{\sin\left( \overset{\_}{{\Delta\phi}_{S}} \right)}}{{R_{F}{\cos\left( \overset{\_}{{\Delta\phi}_{S}} \right)}} + R_{B}^{\prime}} \right\rbrack}} & (217)\end{matrix}$

If R_(B)′ turns out to be too large, i.e. if δ_(max)R_(B)>R_(B)′ thenR_(B)′ is re-defined asR_(B)′≡δ_(max)R_(B)   (218)and _(s)′ is calculated from equation (217).

The proposed location procedure incorporating at least some of thelocation methods utilizing CI,TA or and RX information is shown in FIG.17. The procedure can incorporate location methods presented above (CI,CI+TA and CI+RX algorithms). The process of which is described asfollows: Firstly, in step S1 the timing advance and received signallevels measured by the mobile station from the serving and neighbourcells are collected.

Next, the relevant algorithm input parameters and radio networkparameters such as the base transceiver station coordinates for themeasured cells are collected in step S2.

In step S3, the measurements and the network data are analysed to selectwhich measurements have to be used in the location calculation. Cellsthat do not contribute significantly to improving the location accuracyare removed from the measurement set. As a result of the cell selectionprocedure, serving cell information can potentially be removed and/orneighbour cell information can possibly be totally or partially removed.

In step S4, a cell identity/cell identity plus timing advance locationestimate is made. If serving cell information is available, cellidentity or cell identity and timing advance based location algorithmssuch as described hereinbefore can be used to determine a locationestimate and/or a confidence region for it. The resulting locationestimate is called in this document, the “CI/CI+TA Location Estimate”.

In step S5, an RX location estimate is made. If the neighbour cellinformation is available, a CI+RX based location algorithm is applied.This may be any one of algorithms A to F described hereinbefore orindeed any other suitable algorithm. This is used to estimate the mobilestation location. In selecting the CI+RX based algorithms, the followingcriteria may be taken into account:

-   -   cell sectorisation. For omnidirectional cells, CI+RX based        algorithms can be used. For example, algorithms D (equation        154), E (equation 1-56) or F (equation 158) may be used.

For sector cell cells, CI+RX based algorithms such as algorithms A, Band C described hereinbefore can be used directly. Alternatively,algorithms D, E and F can be extended to the case of sector cells asdescribed previously.

-   -   Absolute/relative level observations

Algorithm C described hereinbefore makes use of relative levelobservations, that is the difference in received signal levelmeasurements from pairs of base transceiver stations. It should beappreciated that algorithms A, B, D, E and F use absolute levelobservations.

-   -   Closed form/iterative algorithms

Algorithms A, B, C, D and E are iterative algorithms. Algorithm F is aclosed form algorithm. It determines the mobile station coordinates as aweighted average of the coordinates of the base transceiver stationsfrom which the received signal level measurements have been collected.One analytical expression (equations 163 and 164) and three empiricalapproximations for the weights are also provided as discussed above

-   -   RX-based domain determination

The RX-based algorithms may require determination of the “domain” inwhich the RX-based location estimated is constrained to lie. Such adomain can be defined using, for example, any of the definitionsprovided hereinbefore.

The obtained location estimate is called for convenience “RX locationestimate”.

In step S6 an RX direction location estimate is obtained. If servingcell information and the “RX location estimate” are available, thedirection from the serving base transceiver station to the “RX locationestimate” is calculated. This direction is used, instead of the actualantenna direction, as an input for a CI+TA based algorithm, as explainedabove. The resulting location estimate is called, for convenience, “RXdirection location estimate”.

In step S7, a virtual BTS location estimate is calculated.

If serving cell information and “RX-direction location estimate” areavailable, the coordinates of the “RX-direction location estimate” isused as an additional (omni-directional) neighbour cell for the RX basedalgorithm. This additional neighbour cell and the RX level measurementassociated with it can be called the “virtual base transceiver station”and “virtual RX level” respectively. The value for the virtual RX levelmeasurement can be selected for example as the maximum RX level from themeasured neighbours or as the average RX level from the measuredneighbours. The resulting location estimate is called “the virtual BTSlocation estimate”.

It should be appreciated that a location estimate obtained by any of thelocation methods discussed can be used as an additional or virtual basestation measurement. The virtual base station measurement is associatedwith a virtual measurement, for example the virtual RX level, and theset of real measurements and the virtual measurement is reprocessedusing any of the location methods described.

In step S8, any one of the obtained location estimates: CI/CI+TAlocation estimate; RX location estimate; RX-direction location estimate;RX location estimate; RX-direction location estimate; or virtual BTSlocation estimate can be selected as the location estimate delivered bythe location process.

Embodiments of the present invention are particularly flexible. Forexample, if the TA measurements and/or serving cell information are notavailable, the process is still able to provide locationinformation—such as the RX location estimate. Likewise, if neighbourinformation and/or RX level measurements are not available, the processis still capable of delivering location estimates using the CI/CI+TAlocation estimate.

Embodiments of the present invention enable both the RX level and TAmeasurements to be used effectively. This information is currentlyprovided in the standard GSM networks. The combined use of timingadvance and RX level measurements “the virtual BTS location estimate”gives the most accurate location estimates in the regions of the highestdensity of the mobile users. In embodiments of the present invention,one of the estimates is selected. Where more than one estimate can beobtained, an averaging or alternatively a weighting of the results maybe provided.

Where a number of different location estimates can be determined, onlyone of the methods may be used. Any suitable criteria can be used fordetermining which of the methods is to be used.

Whilst embodiments of the present invention have been describedparticularly in the context of an arrangement where a number ofdifferent location methods are used, it should be appreciated that insome embodiments different ones of the location algorithms describedhereinbefore can be used alone.

1. A method of estimating the location of a mobile device, comprising the steps of: collecting location information; selecting at least one of a plurality of different location methods to provide a location estimate said methods comprising using cell identity information; and providing a location estimate based on the at least one selected location method.
 2. A method as claimed in claim 1 wherein said at least one location method comprises at least one of the following methods: a method using cell identity information; a method using cell identity information and received signal strength; a method using cell identity information and timing advance information; and a method using cell identity information, received signal strength information and timing advance information.
 3. A method as claimed in claim 1, comprising the step of determining a virtual base station estimate.
 4. A method as claimed in claim 2, further comprising the step of determining a virtual base station estimate, wherein said virtual base station estimate is determined using at least one of the methods of claim
 2. 5. A method as claimed in claim 3, wherein said virtual base station location estimate coupled with at least one virtual measurement and at least one real measurement and said at least one virtual measurement is processed using a location method.
 6. A method as claimed in claim 2, wherein said virtual base station location estimate coupled with at least one virtual measurement and at least one real measurement and said at least one virtual measurement is processed using a location method, and wherein the at least one real and the at least one virtual measurements are processed using a location method as defined in claim
 2. 7. A method as claimed in claim 5, wherein a value for the virtual measurement is one of measured levels, a combination of measured levels, and an average of measured levels.
 8. A method as claimed in claim 1, wherein said at least one location method is selected in dependence on the location information available.
 9. A method as claimed in claim 1, wherein a plurality of location estimates are determined and at least one estimate is used to provide said location estimate.
 10. A method as claimed in claim 1, wherein said location information is collected by said mobile device.
 11. A method as claimed in claim 10, wherein said mobile device is arranged to measure a level of at least one type of information.
 12. A method as claimed in claim 1, wherein said location information comprises at least one of timing advance information and received signal level.
 13. A method as claimed in claim 12, wherein said received signal level is an absolute received signal level or relative received signal level.
 14. A method as claimed in claim 1, wherein said mobile device is in a cellular communications device.
 15. A method as claimed in claim 14, wherein said information is collected for a serving cell of the mobile device.
 16. A method as claimed in claim 14, wherein said information is collected for at least one neighbouring cell.
 17. A method as claimed in claim 14, comprising the step of selecting the or each cell in respect of which location information is collected.
 18. A method as claimed in claim 1, wherein a location estimate is provided using the following algorithm Calculate the total attenuation experienced by a signal transmitted by the i-th BTS while propagating toward a mobile station where i-th level observation is L^(i)) by subtracting from the i-th measured received power, P_(r) ^(i), the maximum power radiated by the i-th BTS, P_(t,max) ^(i): L ^(i) =P _(r) ^(i) −P _(t,max) ^(i) ; i=1, . . . ,N   (11) Stack the level observations from N BTS's in vector L: L=[L¹, . . . ,L^(N)]^(T)   (12) Solve the minimization problem: $\begin{matrix} {\begin{bmatrix} {\hat{\sigma}}_{u}^{2} \\ \hat{x} \\ \hat{y} \end{bmatrix} = {\arg\quad{\min\limits_{\lbrack\begin{matrix} \sigma_{u}^{2} \\ x \\ y \end{matrix}\rbrack}{F\left( {x,{y;\sigma_{u}^{2}}} \right)}}}} & (13) \end{matrix}$ where the cost function F (x,y; σ_(u) ²) is defined as follows: $\begin{matrix} {{{F\left( {x,{y;\sigma_{u}^{2}}} \right)} = {{\ln\quad\sigma_{u}^{2}} + {\ln{{r_{L}\left( {x,y} \right)}}} + {{\frac{1}{\sigma_{u}^{2}}\left\lbrack {L - {m_{L}\left( {x,y} \right)}} \right\rbrack}^{T}{{r_{L}^{- 1}\left( {x,y} \right)}\left\lbrack {L - {m_{L}\left( {x,y} \right)}} \right\rbrack}}}}{and}} & (14) \\ {{m_{L}\left( {x,y} \right)} = \left\lbrack {{\mu_{L}^{1}\left( {x,y} \right)},\ldots\quad,{\mu_{L}^{N}\left( {x,y} \right)}} \right\rbrack^{T}} & (15) \\ {{\mu_{L}^{i}\left( {x,y} \right)} = {{- {{PL}^{i}\left( {d^{i}\left( {x,y} \right)} \right)}} - {{AP}_{tr}^{i}\left( {\psi^{i}\left( {x,y} \right)} \right)}}} & (16) \\ {\left\lbrack {r_{L}\left( {x,y} \right)} \right\rbrack_{ij} = \left\{ {{\begin{matrix} 1 & {i = j} \\ {\rho_{u}^{i,j}\left( {x,y} \right)} & {i \neq j} \end{matrix}\quad i},{j = 1},\ldots\quad,N} \right.} & (17) \end{matrix}$
 19. A method as claimed in claim 1, wherein a location estimate is provided using the following algorithm Calculate the total attenuation experienced by a signal transmitted by the i-th BTS while propagating toward a mobile station where the i-th level observation is L^(i) by subtracting from the i-th measured received power, P_(r) ^(i), the maximum power radiated by the i-th BTS, P_(t,max) ^(i): L ^(i) =P _(r) ^(i) −P _(t,max) ^(i) ; i=1, . . . ,N   (18) Stack level observations from N BTS's in vector L: L=[L¹, . . . ,L^(N)]^(T)   (19) Solve the minimization problem: $\begin{matrix} {\begin{bmatrix} \hat{x} \\ \hat{y} \end{bmatrix} = {\arg\quad{\min\limits_{{\lbrack\begin{matrix} x \\ y \end{matrix}\rbrack}\varepsilon\quad D_{xy}}{F\left( {x,y} \right)}}}} & (20) \end{matrix}$ where the cost function F(x,y) is defined as follows: $\begin{matrix} {{F\left( {x,y} \right)} = {\sum\limits_{i = 1}^{N}\left( {L^{i} + {{PL}^{i}\left( {x,y} \right)} + {{AP}_{tr}^{i}\left( {x,y} \right)}} \right)^{2}}} & (21) \end{matrix}$ and D_(xy) is the domain of existence of x and y. Calculate {circumflex over (σ)}_(u) ² as {circumflex over (σ)}_(u) ² =F({circumflex over (x)},ŷ)   (22)
 20. Amended) A method as claimed in claim 1, wherein a location estimate is provided using the following algorithm: Calculate the total attenuation experienced by a signal transmitted by the i-th BTS while propagating toward a mobile station where the i-th level observation is L^(i)) by subtracting from the i-th measured received power, Pt, the maximum power radiated by the i-th BTS, P_(t,max) ^(i): L ^(i) =P _(r) ^(i) −P _(t,max) ^(i) ; i=1, . . . ,N   (23) Calculate the j-th level difference observation by subtracting the j-th level observation from the level observation L¹ taken as reference: D ^(i) =L ¹ −L ^(j) ; j=2, . . . ,N   (24) Stack the N−1 difference of level observations in a vector D: D=[D², . . . ,D^(N)]^(T)   (25) Solve the minimization problem $\begin{matrix} {{\begin{bmatrix} \hat{x} \\ \hat{y} \end{bmatrix} = {\arg\quad{\min\limits_{{\lbrack\begin{matrix} x \\ y \end{matrix}\rbrack}\varepsilon\quad D_{xy}}{F\left( {x,y} \right)}}}}{where}} & (26) \\ {{{F\left( {x,y} \right)} = {{\sum\limits_{j = 2}^{N}\left( {D^{j} - {\mu_{D}^{j}\left( {x,y} \right)}} \right)^{2}} - {\frac{1}{N}\left( {{\sum\limits_{j = 2}^{N}D^{j}} - {\mu_{D}^{j}\left( {x,y} \right)}} \right)^{2}}}}{and}} & (27) \\ {{\mu_{D}^{j}\left( {x,y} \right)} = {{- \left\lbrack {{{PL}^{1}\left( {d^{1}\left( {x,y} \right)} \right)} - {{PL}^{j}\left( {d^{j}\left( {x,y} \right)} \right)}} \right\rbrack} - \left\lbrack {{{AP}_{tr}^{1}\left( {\psi^{1}\left( {x,y} \right)} \right)} - {{AP}_{tr}^{j}\left( {\psi^{j}\left( {x,y} \right)} \right)}} \right\rbrack}} & (28) \end{matrix}$ D_(xy) is the domain of existence of x and y.
 21. A method as claimed in claim 1, wherein a location estimate is provided using an algorithm solving the following equation in x and y: $\left\{ {\begin{matrix} {{\sum\limits_{i = 1}^{N}{{F^{i}\left( {x,y} \right)}\left( {x - x^{i}} \right)}} = 0} \\ {{\sum\limits_{i = 1}^{N}{{F^{i}\left( {x,y} \right)}\left( {y - y^{i}} \right)}} = 0} \end{matrix};{{\left( {x,y} \right) \in {{??}\quad{where}{F^{i}\left( {x,y} \right)}}} = {\frac{2{B^{i}/{C^{i}\left( d_{0} \right)}}}{\left( {2\pi} \right)^{3/2}\sigma_{u}^{i}\ln\quad 10}\quad{\frac{\exp{\left\{ {{- \frac{1}{2\sigma_{u}^{i^{2}}}}\left( {{B^{i}\log_{10}{d^{i}\left( {x,y} \right)}} - z^{i} + A^{i}} \right)^{2}} \right\} \cdot}}{\left\lbrack {d^{i}\left( {x,y} \right)} \right\rbrack^{4}}\quad\left\lbrack {\frac{B^{i}\left( {{B^{i}\log_{10}{d^{i}\left( {x,y} \right)}} - x^{i} + A^{i}} \right)}{2\sigma_{u}^{i^{2}}\ln\quad 10} - 1} \right\rbrack}}}} \right.$
 22. A method as claimed in claim 1, wherein a location estimate is provided using an algorithm solving the following equation in x and y: $\left\{ {\begin{matrix} {{\sum\limits_{i = 1}^{N}\left\lbrack {{{- \frac{\mathcal{I}_{i}}{\Re }}\left( {x - x^{i}} \right)} - {\frac{\left( {{\overset{\sim}{\mathcal{I}}}_{i} - 1} \right)}{\Re }\left\{ {{\left( x^{i} \right)^{2}x} - {x^{i}{y^{i}\left( {y - y^{i}} \right)}}} \right\}}} \right\rbrack} = 0} \\ {{\sum\limits_{i = 1}^{N}\left\lbrack {{{- \frac{\mathcal{I}_{i}}{\Re }}\left( {y - y^{i}} \right)} - {\frac{\left( {{\overset{\sim}{\mathcal{I}}}_{i} - 1} \right)}{\Re }\left\{ {{\left( y^{i} \right)^{2}y} - {x^{i}{y^{i}\left( {x - x^{i}} \right)}}} \right\}}} \right\rbrack} = 0} \end{matrix};{\left( {x,y} \right) \in {??}}} \right.$
 23. A method as claimed in claim 1, wherein a location estimate is provided using an algorithm based on the following equation: ${\hat{x} = \frac{\sum\limits_{i = 1}^{N}\frac{x^{i}}{\mathcal{I}_{i0}}}{\sum\limits_{i = 1}^{N}\frac{1}{\mathcal{I}_{i0}}}};{\hat{y} = \frac{\sum\limits_{i = 1}^{N}\frac{y^{i}}{\mathcal{I}_{i0}}}{\sum\limits_{i = 1}^{N}\frac{1}{\mathcal{I}_{i0}}}};{\left( {\hat{x},\hat{y}} \right) \in {??}}$
 24. A method as claimed in claim 1, wherein said location estimate is provided by one of a iterative and a closed form method.
 25. A method as claimed in claim 1, wherein said location estimate is provided by one of a linear and non linear method.
 26. A system for estimating the location of a mobile device, comprising: means for collecting location information; means for selecting at least one of a plurality of different location methods to provide a location estimate said methods using cell identity information; and means for providing a location estimate based on the at least one selected location method. 